Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+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
parallel_in_time.h
Go to the documentation of this file.
1
317 *  
318 *   #include "BraidFuncs.hh"
319 *  
320 * @endcode
321 *
322 * This advances the solution forward by one time step.
323 * First some data is collected from the status struct,
324 * namely the start and stop time and the current timestep
325 * number. The timestep size @f$\Delta t@f$ is calculated,
326 * and the step function from the HeatEquation is used to
327 * advance the solution.
328 *
329 * @code
330 *   int my_Step(braid_App app,
331 *   braid_Vector ustop,
332 *   braid_Vector fstop,
333 *   braid_Vector u,
334 *   braid_StepStatus status)
335 *   {
336 *   UNUSED(ustop);
337 *   UNUSED(fstop);
338 *   double tstart; /* current time */
339 *   double tstop; /* evolve to this time*/
340 *   int level;
341 *   double deltaT;
342 *  
343 *   int index;
344 *   braid_StepStatusGetLevel(status, &level);
345 *   braid_StepStatusGetTstartTstop(status, &tstart, &tstop);
346 *   braid_StepStatusGetTIndex(status, &index);
347 *  
348 *   deltaT = tstop - tstart;
349 *  
350 *   ::Vector<double>& solution = u->data;
351 *  
352 *   HeatEquation<2>& heateq = app->eq;
353 *  
354 *   heateq.step(solution, deltaT, tstart, index);
355 *  
356 *   return 0;
357 *   }
358 *  
359 *  
360 * @endcode
361 *
362 * In this function we initialize a vector at an arbitrary time.
363 * At this point we don't know anything about what the solution
364 * looks like, and we can really initialize to anything, so in
365 * this case use reinit to initialize the memory and set the
366 * values to zero.
367 *
368 * @code
369 *   int
370 *   my_Init(braid_App app,
371 *   double t,
372 *   braid_Vector *u_ptr)
373 *   {
374 *   my_Vector *u = new(my_Vector);
375 *   int size = app->eq.size();
376 *   u->data.reinit(size);
377 *  
378 *   app->eq.initialize(t, u->data);
379 *  
380 *   *u_ptr = u;
381 *  
382 *   return 0;
383 *   }
384 *  
385 * @endcode
386 *
387 * Here we need to copy the vector u into the vector v. We do this
388 * by allocating a new vector, then reinitializing the deal.ii
389 * vector to the correct size. The deal.ii reinitialization sets
390 * every value to zero, so next we need to iterate over the vector
391 * u and copy the values to the new vector v.
392 *
393 * @code
394 *   int
395 *   my_Clone(braid_App app,
396 *   braid_Vector u,
397 *   braid_Vector *v_ptr)
398 *   {
399 *   UNUSED(app);
400 *   my_Vector *v = new(my_Vector);
401 *   int size = u->data.size();
402 *   v->data.reinit(size);
403 *   for(size_t i=0, end=v->data.size(); i != end; ++i)
404 *   {
405 *   v->data[i] = u->data[i];
406 *   }
407 *   *v_ptr = v;
408 *  
409 *   return 0;
410 *   }
411 *  
412 * @endcode
413 *
414 * Here we need to free the memory used by vector u. This is
415 * pretty simple since the deal.ii vector is stored inside the
416 * XBraid vector, so we just delete the XBraid vector u and it
417 * puts the deal.ii vector out of scope and releases its memory.
418 *
419 * @code
420 *   int
421 *   my_Free(braid_App app,
422 *   braid_Vector u)
423 *   {
424 *   UNUSED(app);
425 *   delete u;
426 *  
427 *   return 0;
428 *   }
429 *  
430 * @endcode
431 *
432 * This is to perform an axpy type operation. That is to say we
433 * do @f$y = \alpha x + \beta y@f$. Fortunately deal.ii already has
434 * this operation built in to its vector class, so we get the
435 * reference to the vector y and call the sadd method.
436 *
437 * @code
438 *   int my_Sum(braid_App app,
439 *   double alpha,
440 *   braid_Vector x,
441 *   double beta,
442 *   braid_Vector y)
443 *   {
444 *   UNUSED(app);
445 *   Vector<double>& vec = y->data;
446 *   vec.sadd(beta, alpha, x->data);
447 *  
448 *   return 0;
449 *   }
450 *  
451 * @endcode
452 *
453 * This calculates the spatial norm using the l2 norm. According
454 * to XBraid, this could be just about any spatial norm but we'll
455 * keep it simple and used deal.ii vector's built in l2_norm method.
456 *
457 * @code
458 *   int
459 *   my_SpatialNorm(braid_App app,
460 *   braid_Vector u,
461 *   double *norm_ptr)
462 *   {
463 *   UNUSED(app);
464 *   double dot = 0.0;
465 *   dot = u->data.l2_norm();
466 *   *norm_ptr = dot;
467 *  
468 *   return 0;
469 *   }
470 *  
471 * @endcode
472 *
473 * This function is called at various points depending on the access
474 * level specified when configuring the XBraid struct. This function
475 * is used to print out data during the run time, such as plots of the
476 * data. The status struct contains a ton of information about the
477 * simulation run. Here we get the current time and timestep number.
478 * The output_results function is called to plot the solution data.
479 * If the method of manufactured solutions is being used, then the
480 * error of this time step is computed and processed.
481 *
482 * @code
483 *   int
484 *   my_Access(braid_App app,
485 *   braid_Vector u,
486 *   braid_AccessStatus astatus)
487 *   {
488 *   double t;
489 *   int index;
490 *  
491 *   braid_AccessStatusGetT(astatus, &t);
492 *   braid_AccessStatusGetTIndex(astatus, &index);
493 *  
494 *   app->eq.output_results(index, t, u->data);
495 *  
496 *   #if DO_MFG
497 *   if(index == app->final_step)
498 *   {
499 *   app->eq.process_solution(t, index, u->data);
500 *   }
501 *   #endif
502 *  
503 *   return 0;
504 *   }
505 *  
506 * @endcode
507 *
508 * This calculates the size of buffer needed to pack the solution
509 * data into a linear buffer for transfer to another processor via
510 * MPI. We query the size of the data from the HeatEquation class
511 * and return the buffer size.
512 *
513 * @code
514 *   int
515 *   my_BufSize(braid_App app,
516 *   int *size_ptr,
517 *   braid_BufferStatus bstatus)
518 *   {
519 *   UNUSED(bstatus);
520 *   int size = app->eq.size();
521 *   *size_ptr = (size+1)*sizeof(double);
522 *  
523 *   return 0;
524 *   }
525 *  
526 * @endcode
527 *
528 * This function packs a linear buffer with data so that the buffer
529 * may be sent to another processor via MPI. The buffer is cast to
530 * a type we can work with. The first element of the buffer is the
531 * size of the buffer. Then we iterate over soltuion vector u and
532 * fill the buffer with our solution data. Finally we tell XBraid
533 * how much data we wrote.
534 *
535 * @code
536 *   int
537 *   my_BufPack(braid_App app,
538 *   braid_Vector u,
539 *   void *buffer,
540 *   braid_BufferStatus bstatus)
541 *   {
542 *  
543 *   UNUSED(app);
544 *   double *dbuffer = (double*)buffer;
545 *   int size = u->data.size();
546 *   dbuffer[0] = size;
547 *   for(int i=0; i != size; ++i)
548 *   {
549 *   dbuffer[i+1] = (u->data)[i];
550 *   }
551 *   braid_BufferStatusSetSize(bstatus, (size+1)*sizeof(double));
552 *  
553 *   return 0;
554 *   }
555 *  
556 * @endcode
557 *
558 * This function unpacks a buffer that was recieved from a different
559 * processor via MPI. The size of the buffer is read from the first
560 * element, then we iterate over the size of the buffer and fill
561 * the values of solution vector u with the data in the buffer.
562 *
563 * @code
564 *   int
565 *   my_BufUnpack(braid_App app,
566 *   void *buffer,
567 *   braid_Vector *u_ptr,
568 *   braid_BufferStatus bstatus)
569 *   {
570 *   UNUSED(app);
571 *   UNUSED(bstatus);
572 *  
573 *   my_Vector *u = NULL;
574 *   double *dbuffer = (double*)buffer;
575 *   int size = static_cast<int>(dbuffer[0]);
576 *   u = new(my_Vector);
577 *   u->data.reinit(size);
578 *  
579 *   for(int i = 0; i != size; ++i)
580 *   {
581 *   (u->data)[i] = dbuffer[i+1];
582 *   }
583 *   *u_ptr = u;
584 *  
585 *   return 0;
586 *   }
587 * @endcode
588
589
590<a name="ann-src/BraidFuncs.hh"></a>
591<h1>Annotated version of src/BraidFuncs.hh</h1>
592 *
593 *
594 *
595 *
596 * @code
597 *   /* -----------------------------------------------------------------------------
598 *   *
599 *   * SPDX-License-Identifier: LGPL-2.1-or-later
600 *   * Copyright (C) 2018 by Joshua Christopher
601 *   *
602 *   * This file is part of the deal.II code gallery.
603 *   *
604 *   * -----------------------------------------------------------------------------
605 *   */
606 *  
607 *   #ifndef _BRAIDFUNCS_H_
608 *   #define _BRAIDFUNCS_H_
609 *  
610 *   /**
611 *   * \file BraidFuncs.cc
612 *   * \brief Contains the implementation of the mandatory X-Braid functions
613 *   *
614 *   * X-Braid mandates several functions in order to drive the solution.
615 *   * This file contains the implementation of said mandatory functions.
616 *   * See the X-Braid documentation for more information.
617 *   * There are several functions that are optional in X-Braid that may
618 *   * or may not be implemented in here.
619 *   *
620 *   */
621 *  
622 *  
623 *   /*-------- Third Party --------*/
624 *   #include <deal.II/numerics/vector_tools.h>
625 *  
626 *   #include <braid.h>
627 *   #include <braid_test.h>
628 *  
629 *   /*-------- Project --------*/
630 *   #include "HeatEquation.hh"
631 *  
632 * @endcode
633 *
634 * This struct contains all data that changes with time. For now
635 * this is just the solution data. When doing AMR this should
636 * probably include the triangulization, the sparsity patter,
637 * constraints, etc.
638 *
639 * @code
640 *   /**
641 *   * \brief Struct that contains the deal.ii vector.
642 *   */
643 *   typedef struct _braid_Vector_struct
644 *   {
645 *   ::Vector<double> data;
646 *   } my_Vector;
647 *  
648 * @endcode
649 *
650 * This struct contains all the data that is unchanging with time.
651 *
652 * @code
653 *   /**
654 *   * \brief Struct that contains the HeatEquation and final
655 *   * time step number.
656 *   */
657 *   typedef struct _braid_App_struct
658 *   {
659 *   HeatEquation<2> eq;
660 *   int final_step;
661 *   } my_App;
662 *  
663 *  
664 *   /**
665 *   * @brief my_Step - Takes a step in time, advancing the u vector
666 *   *
667 *   * @param app - The braid app struct
668 *   * @param ustop - The solution data at the end of this time step
669 *   * @param fstop - RHS data (such as forcing function?)
670 *   * @param u - The solution data at the beginning of this time step
671 *   * @param status - Status structure that contains various info of this time
672 *   *
673 *   * @return Success (0) or failure (1)
674 *   **/
675 *   int my_Step(braid_App app,
676 *   braid_Vector ustop,
677 *   braid_Vector fstop,
678 *   braid_Vector u,
679 *   braid_StepStatus status);
680 *  
681 *  
682 *   /**
683 *   * @brief my_Init - Initializes a solution data at the given time
684 *   * For now, initializes the solution to zero no matter what time we are at
685 *   *
686 *   * @param app - The braid app struct containing user data
687 *   * @param t - Time at which the solution is initialized
688 *   * @param u_ptr - The solution data that needs to be filled
689 *   *
690 *   * @return Success (0) or failure (1)
691 *   **/
692 *   int
693 *   my_Init(braid_App app,
694 *   double t,
695 *   braid_Vector *u_ptr);
696 *  
697 *  
698 *   /**
699 *   * @brief my_Clone - Clones a vector into a new vector
700 *   *
701 *   * @param app - The braid app struct containing user data
702 *   * @param u - The existing vector containing data
703 *   * @param v_ptr - The empty vector that needs to be filled
704 *   *
705 *   * @return Success (0) or failure (1)
706 *   **/
707 *   int
708 *   my_Clone(braid_App app,
709 *   braid_Vector u,
710 *   braid_Vector *v_ptr);
711 *  
712 *  
713 *   /**
714 *   * @brief my_Free - Deletes a vector
715 *   *
716 *   * @param app - The braid app struct containing user data
717 *   * @param u - The vector that needs to be deleted
718 *   *
719 *   * @return Success (0) or failure (1)
720 *   **/
721 *   int
722 *   my_Free(braid_App app,
723 *   braid_Vector u);
724 *  
725 *  
726 *   /**
727 *   * @brief my_Sum - Sums two vectors in an AXPY operation
728 *   * The operation is y = alpha*x + beta*y
729 *   *
730 *   * @param app - The braid app struct containing user data
731 *   * @param alpha - The coefficient in front of x
732 *   * @param x - A vector that is multiplied by alpha then added to y
733 *   * @param beta - The coefficient of y
734 *   * @param y - A vector that is multiplied by beta then summed with x
735 *   *
736 *   * @return Success (0) or failure (1)
737 *   **/
738 *   int
739 *   my_Sum(braid_App app,
740 *   double alpha,
741 *   braid_Vector x,
742 *   double beta,
743 *   braid_Vector y);
744 *  
745 *   /**
746 *   * \brief Returns the spatial norm of the provided vector
747 *   *
748 *   * Calculates and returns the spatial norm of the provided vector.
749 *   * Interestingly enough, X-Braid does not specify a particular norm.
750 *   * to keep things simple, we implement the Euclidean norm.
751 *   *
752 *   * \param app - The braid app struct containing user data
753 *   * \param u - The vector we need to take the norm of
754 *   * \param norm_ptr - Pointer to the norm that was calculated, need to modify this
755 *   * \return Success (0) or failure (1)
756 *   */
757 *   int
758 *   my_SpatialNorm(braid_App app,
759 *   braid_Vector u,
760 *   double *norm_ptr);
761 *  
762 *   /**
763 *   * \brief Allows the user to output details
764 *   *
765 *   * The Access function is called at various points to allow the user to output
766 *   * information to the screen or to files.
767 *   * The astatus parameter provides various information about the simulation,
768 *   * see the XBraid documentation for details on what information you can get.
769 *   * Example information is what the current timestep number and current time is.
770 *   * If the access level (in parallel_in_time.cc) is set to 0, this function is
771 *   * never called.
772 *   * If the access level is set to 1, the function is called after the last
773 *   * XBraid cycle.
774 *   * If the access level is set to 2, it is called every XBraid cycle.
775 *   *
776 *   * \param app - The braid app struct containing user data
777 *   * \param u - The vector containing the data at the status provided
778 *   * \param astatus - The Braid status structure
779 *   * \return Success (0) or failure (1)
780 *   */
781 *   int
782 *   my_Access(braid_App app,
783 *   braid_Vector u,
784 *   braid_AccessStatus astatus);
785 *  
786 *   /**
787 *   * \brief Calculates the size of a buffer for MPI data transfer
788 *   *
789 *   * Calculates the size of the buffer that is needed to transfer
790 *   * a solution vector to another processor.
791 *   * The bstatus parameter provides various information on the
792 *   * simulation, see the XBraid documentation for all possible
793 *   * fields.
794 *   *
795 *   * \param app - The braid app struct containing user data
796 *   * \param size_ptr A pointer to the calculated size
797 *   * \param bstatus The XBraid status structure
798 *   * \return Success (0) or failure (1)
799 *   */
800 *   int
801 *   my_BufSize(braid_App app,
802 *   int *size_ptr,
803 *   braid_BufferStatus bstatus);
804 *  
805 *   /**
806 *   * \brief Linearizes a vector to be sent to another processor
807 *   *
808 *   * Linearizes (packs) a data buffer with the contents of
809 *   * some solution state u.
810 *   *
811 *   * \param app - The braid app struct containing user data
812 *   * \param u The vector that must be packed into buffer
813 *   * \param buffer The buffer that must be filled with u
814 *   * \param bstatus The XBraid status structure
815 *   * \return Success (0) or failure (1)
816 *   */
817 *   int
818 *   my_BufPack(braid_App app,
819 *   braid_Vector u,
820 *   void *buffer,
821 *   braid_BufferStatus bstatus);
822 *  
823 *   /**
824 *   * \brief Unpacks a vector that was sent from another processor
825 *   *
826 *   * Unpacks a linear data buffer into the vector pointed to by
827 *   * u_ptr.
828 *   *
829 *   * \param app - The braid app struct containing user data
830 *   * \param buffer The buffer that must be unpacked
831 *   * \param u_ptr The pointer to the vector that is filled
832 *   * \param bstatus The XBraid status structure
833 *   * \return Success (0) or failure (1)
834 *   */
835 *   int
836 *   my_BufUnpack(braid_App app,
837 *   void *buffer,
838 *   braid_Vector *u_ptr,
839 *   braid_BufferStatus bstatus);
840 *  
841 *   #endif // _BRAIDFUNCS_H_
842 * @endcode
843
844
845<a name="ann-src/HeatEquation.hh"></a>
846<h1>Annotated version of src/HeatEquation.hh</h1>
847 *
848 *
849 *
850 *
851 * @code
852 *   /* -----------------------------------------------------------------------------
853 *   *
854 *   * SPDX-License-Identifier: LGPL-2.1-or-later
855 *   * Copyright (C) 2018 by Joshua Christopher
856 *   *
857 *   * This file is part of the deal.II code gallery.
858 *   *
859 *   * -----------------------------------------------------------------------------
860 *   */
861 *  
862 *   #ifndef _HEATEQUATION_H_
863 *   #define _HEATEQUATION_H_
864 *  
865 *   #include <deal.II/base/utilities.h>
866 *   #include <deal.II/base/quadrature_lib.h>
867 *   #include <deal.II/base/function.h>
868 *   #include <deal.II/base/logstream.h>
869 *   #include <deal.II/lac/vector.h>
870 *   #include <deal.II/lac/full_matrix.h>
871 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
872 *   #include <deal.II/lac/sparse_matrix.h>
873 *   #include <deal.II/lac/solver_cg.h>
874 *   #include <deal.II/lac/precondition.h>
875 *   #include <deal.II/lac/affine_constraints.h>
876 *   #include <deal.II/grid/tria.h>
877 *   #include <deal.II/grid/grid_generator.h>
878 *   #include <deal.II/grid/grid_refinement.h>
879 *   #include <deal.II/grid/grid_out.h>
880 *   #include <deal.II/grid/tria_accessor.h>
881 *   #include <deal.II/grid/tria_iterator.h>
882 *   #include <deal.II/dofs/dof_handler.h>
883 *   #include <deal.II/dofs/dof_accessor.h>
884 *   #include <deal.II/dofs/dof_tools.h>
885 *   #include <deal.II/fe/fe_q.h>
886 *   #include <deal.II/fe/fe_values.h>
887 *   #include <deal.II/numerics/data_out.h>
888 *   #include <deal.II/numerics/vector_tools.h>
889 *   #include <deal.II/numerics/error_estimator.h>
890 *   #include <deal.II/numerics/solution_transfer.h>
891 *   #include <deal.II/numerics/matrix_tools.h>
892 *   #include <deal.II/base/convergence_table.h>
893 *  
894 *   #include <fstream>
895 *  
896 *   using namespace dealii;
897 *  
898 * @endcode
899 *
900 * The HeatEquation class is describes the finite element
901 * solver for the heat equation. It contains all the functions
902 * needed to define the problem domain and advance the solution
903 * in time.
904 *
905 * @code
906 *   template <int dim>
907 *   class HeatEquation
908 *   {
909 *   public:
910 *   HeatEquation();
911 *   void define();
912 *   void step(Vector<double>& braid_data,
913 *   double deltaT,
914 *   double a_time,
915 *   int a_time_idx);
916 *  
917 *   int size() const; /// Returns the size of the solution vector
918 *  
919 *   void output_results(int a_time_idx,
920 *   double a_time,
921 *   Vector<double>& a_solution) const;
922 *  
923 *   void initialize(double a_time,
924 *   Vector<double>& a_vector) const;
925 *  
926 *   void process_solution(double a_time,
927 *   int a_index,
928 *   const Vector<double>& a_vector);
929 *  
930 *   private:
931 *   void setup_system();
932 *   void solve_time_step(Vector<double>& a_solution);
933 *  
934 *   Triangulation<dim> triangulation;
935 *   FE_Q<dim> fe;
936 *   DoFHandler<dim> dof_handler;
937 *  
938 *   AffineConstraints<double> constraints;
939 *  
940 *   SparsityPattern sparsity_pattern;
941 *   SparseMatrix<double> mass_matrix;
942 *   SparseMatrix<double> laplace_matrix;
943 *   SparseMatrix<double> system_matrix;
944 *  
945 *   Vector<double> system_rhs;
946 *  
947 *   std::ofstream myfile;
948 *  
949 *   const double theta;
950 *  
951 * @endcode
952 *
953 * These were originally in the run() function but because
954 * I am splitting the run() function up into define and step
955 * they need to become member data
956 *
957 * @code
958 *   Vector<double> tmp;
959 *   Vector<double> forcing_terms;
960 *  
961 *   ConvergenceTable convergence_table;
962 *   };
963 *  
964 * @endcode
965 *
966 * The RightHandSide class describes the RHS of the governing
967 * equations. In this case, it is the forcing function.
968 *
969 * @code
970 *   template <int dim>
971 *   class RightHandSide : public Function<dim>
972 *   {
973 *   public:
974 *   RightHandSide ()
975 *   :
976 *   Function<dim>(),
977 *   period (0.2)
978 *   {}
979 *  
980 *   virtual double value (const Point<dim> &p,
981 *   const unsigned int component = 0) const override;
982 *  
983 *   private:
984 *   const double period;
985 *   };
986 *  
987 * @endcode
988 *
989 * The BoundaryValues class describes the boundary conditions
990 * of the governing equations.
991 *
992 * @code
993 *   template <int dim>
994 *   class BoundaryValues : public Function<dim>
995 *   {
996 *   public:
997 *   virtual double value (const Point<dim> &p,
998 *   const unsigned int component = 0) const override;
999 *   };
1000 *  
1001 * @endcode
1002 *
1003 * The RightHandSideMFG class describes the right hand side
1004 * function when doing the method of manufactured solutions.
1005 *
1006 * @code
1007 *   template <int dim>
1008 *   class RightHandSideMFG : public Function<dim>
1009 *   {
1010 *   public:
1011 *   virtual double value (const Point<dim> &p,
1012 *   const unsigned int component = 0) const override;
1013 *   };
1014 *  
1015 * @endcode
1016 *
1017 * The InitialValuesMFG class describes the initial values
1018 * when doing the method of manufactured solutions.
1019 *
1020 * @code
1021 *   template <int dim>
1022 *   class InitialValuesMFG : public Function<dim>
1023 *   {
1024 *   public:
1025 *   virtual double value (const Point<dim> &p,
1026 *   const unsigned int component = 0) const override;
1027 *   };
1028 *  
1029 * @endcode
1030 *
1031 * Provides the exact value for the manufactured solution. This
1032 * is used for the boundary conditions as well.
1033 *
1034 * @code
1035 *   template <int dim>
1036 *   class ExactValuesMFG : public Function<dim>
1037 *   {
1038 *   public:
1039 *   /**
1040 *   * \brief Computes the value at the given point and member data time
1041 *   *
1042 *   * Computes the exact value of the manufactured solution at point p and
1043 *   * the member data time. See the class documentation and the design doc
1044 *   * for details on what the exact solution is.
1045 *   *
1046 *   * \param p The point that the exact solution is computed at
1047 *   * \param component The component of the exact solution (always 0 for now)
1048 *   * \return double The exact value that was computed
1049 *   */
1050 *   virtual double value (const Point<dim> &p,
1051 *   const unsigned int component = 0) const override;
1052 *  
1053 *   /**
1054 *   * \brief Computes the gradient of the exact solution at the given point
1055 *   *
1056 *   * Computes the gradient of the exact/manufactured solution value at
1057 *   * point p and member data time. See the design doc for details on
1058 *   * what the gradient of the exact solution is
1059 *   *
1060 *   * \param p The point that the gradient is calculated at
1061 *   * \param component The component of the system of equations this gradient is for
1062 *   * \return Tensor<1,dim> A rank 1 tensor that contains the gradient
1063 *   * in each spatial dimension
1064 *   */
1065 *   virtual Tensor<1,dim> gradient (const Point<dim> &p,
1066 *   const unsigned int component = 0) const override;
1067 *   };
1068 *  
1069 *  
1070 *   #include "HeatEquationImplem.hh"
1071 *  
1072 *   #endif // _HEATEQUATION_H_
1073 * @endcode
1074
1075
1076<a name="ann-src/HeatEquationImplem.hh"></a>
1077<h1>Annotated version of src/HeatEquationImplem.hh</h1>
1078 *
1079 *
1080 *
1081 *
1082 * @code
1083 *   /* -----------------------------------------------------------------------------
1084 *   *
1085 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1086 *   * Copyright (C) 2018 by Joshua Christopher
1087 *   *
1088 *   * This file is part of the deal.II code gallery.
1089 *   *
1090 *   * -----------------------------------------------------------------------------
1091 *   */
1092 *  
1093 *   #include "Utilities.hh"
1094 *  
1095 *   #include <iomanip>
1096 *   #include <math.h>
1097 *  
1098 * @endcode
1099 *
1100 * Calculates the forcing function for the RightHandSide. See the
1101 * documentation for the math.
1102 *
1103 * @code
1104 *   template <int dim>
1105 *   double RightHandSide<dim>::value (const Point<dim> &p,
1106 *   const unsigned int component) const
1107 *   {
1108 *   (void) component;
1109 *   Assert (component == 0, ExcIndexRange(component, 0, 1));
1110 *   Assert (dim == 2, ExcNotImplemented());
1111 *  
1112 *   double time = this->get_time();
1113 *  
1114 *   if ((p[0] > 0.5) && (p[1] > -0.5))
1115 *   {
1116 *   return std::exp(-0.5*(time-0.125)*(time-0.125)/(0.005));
1117 *   }
1118 *   else if ((p[0] > -0.5) && (p[1] > 0.5))
1119 *   {
1120 *   return std::exp(-0.5*(time-0.375)*(time-0.375)/(0.005));
1121 *   }
1122 *   else
1123 *   {
1124 *   return 0;
1125 *   }
1126 *  
1127 *   return 0; // No forcing function
1128 *   }
1129 *  
1130 * @endcode
1131 *
1132 * Calculates the forcing function for the method of manufactured
1133 * solutions. See the documentation for the math.
1134 *
1135 * @code
1136 *   template <int dim>
1137 *   double RightHandSideMFG<dim>::value (const Point<dim> &p,
1138 *   const unsigned int component) const
1139 *   {
1140 *   (void) component;
1141 *   Assert (component == 0, ExcIndexRange(component, 0, 1));
1142 *   Assert (dim == 2, ExcNotImplemented());
1143 *  
1144 *   double time = this->get_time();
1145 *  
1146 *   double pi = numbers::PI;
1147 *   return 4*pi*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
1148 *   }
1149 *  
1150 * @endcode
1151 *
1152 * Calculates the boundary conditions, essentially zero everywhere.
1153 *
1154 * @code
1155 *   template <int dim>
1156 *   double BoundaryValues<dim>::value (const Point<dim> &p,
1157 *   const unsigned int component) const
1158 *   {
1159 *   UNUSED(p);
1160 *   (void) component;
1161 *   Assert (component == 0, ExcIndexRange(component, 0, 1));
1162 *   return 0;
1163 *   }
1164 *  
1165 * @endcode
1166 *
1167 * Calculates the exact solution (and thus also boundary conditions)
1168 * for the method of manufactured solutions.
1169 *
1170 * @code
1171 *   template <int dim>
1172 *   double ExactValuesMFG<dim>::value (const Point<dim> &p,
1173 *   const unsigned int component) const
1174 *   {
1175 *   (void) component;
1176 *   Assert (component == 0, ExcIndexRange(component, 0, 1));
1177 *  
1178 *   double time = this->get_time();
1179 *   const double pi = numbers::PI;
1180 *  
1181 *   return std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
1182 *   }
1183 *  
1184 * @endcode
1185 *
1186 * Calculates the gradient of the exact solution for the method of manufactured
1187 * solutions. See the documentation for the math.
1188 *
1189 * @code
1190 *   template <int dim>
1191 *   Tensor<1,dim> ExactValuesMFG<dim>::gradient (const Point<dim> &p,
1192 *   const unsigned int) const
1193 *   {
1194 *   Assert (dim == 2, ExcNotImplemented());
1195 *  
1196 *   Tensor<1,dim> return_value;
1197 *   const double pi = numbers::PI;
1198 *   double time = this->get_time();
1199 *   return_value[0] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[1])*std::sin(2*pi*p[0]);
1200 *   return_value[1] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::sin(2*pi*p[1]);
1201 *   return return_value;
1202 *   }
1203 *  
1204 * @endcode
1205 *
1206 * Calculates the initial values for the method of manufactured solutions.
1207 * See the documentation for the math.
1208 *
1209 * @code
1210 *   template <int dim>
1211 *   double InitialValuesMFG<dim>::value (const Point<dim> &p,
1212 *   const unsigned int component) const
1213 *   {
1214 *   (void) component;
1215 *   Assert (component == 0, ExcIndexRange(component, 0, 1));
1216 *   const double pi = numbers::PI;
1217 *  
1218 *   return std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
1219 *   }
1220 *  
1221 *   template <int dim>
1222 *   HeatEquation<dim>::HeatEquation ()
1223 *   :
1224 *   fe(1),
1225 *   dof_handler(triangulation),
1226 *   theta(0.5)
1227 *   {
1228 *   }
1229 *  
1230 *   template <int dim>
1231 *   void HeatEquation<dim>::initialize(double a_time,
1232 *   Vector<double>& a_vector) const
1233 *   {
1234 *   #if DO_MFG
1235 * @endcode
1236 *
1237 * We only initialize values in the manufactured solution case
1238 *
1239 * @code
1240 *   InitialValuesMFG<dim> iv_function;
1241 *   iv_function.set_time(a_time);
1242 *   VectorTools::project (dof_handler, constraints,
1243 *   QGauss<dim>(fe.degree+1), iv_function,
1244 *   a_vector);
1245 *   #else
1246 *   UNUSED(a_time);
1247 *   UNUSED(a_vector);
1248 *   #endif // DO_MFG
1249 * @endcode
1250 *
1251 * If not the MFG solution case, a_vector is already zero'd so do nothing
1252 *
1253 * @code
1254 *   }
1255 *  
1256 *   template <int dim>
1257 *   void HeatEquation<dim>::setup_system()
1258 *   {
1259 *   dof_handler.distribute_dofs(fe);
1260 *  
1261 *   constraints.clear ();
1263 *   constraints);
1264 *   constraints.close();
1265 *  
1266 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
1267 *   DoFTools::make_sparsity_pattern(dof_handler,
1268 *   dsp,
1269 *   constraints,
1270 *   /*keep_constrained_dofs = */ true);
1271 *   sparsity_pattern.copy_from(dsp);
1272 *  
1273 *   mass_matrix.reinit(sparsity_pattern);
1274 *   laplace_matrix.reinit(sparsity_pattern);
1275 *   system_matrix.reinit(sparsity_pattern);
1276 *  
1277 *   MatrixCreator::create_mass_matrix(dof_handler,
1278 *   QGauss<dim>(fe.degree+1),
1279 *   mass_matrix);
1280 *   MatrixCreator::create_laplace_matrix(dof_handler,
1281 *   QGauss<dim>(fe.degree+1),
1282 *   laplace_matrix);
1283 *  
1284 *   system_rhs.reinit(dof_handler.n_dofs());
1285 *   }
1286 *  
1287 *  
1288 *   template <int dim>
1289 *   void HeatEquation<dim>::solve_time_step(Vector<double>& a_solution)
1290 *   {
1291 *   SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
1292 *   SolverCG<> cg(solver_control);
1293 *  
1294 *   PreconditionSSOR<> preconditioner;
1295 *   preconditioner.initialize(system_matrix, 1.0);
1296 *  
1297 *   cg.solve(system_matrix, a_solution, system_rhs,
1298 *   preconditioner);
1299 *  
1300 *   constraints.distribute(a_solution);
1301 *   }
1302 *  
1303 *  
1304 *  
1305 *   template <int dim>
1306 *   void HeatEquation<dim>::output_results(int a_time_idx,
1307 *   double a_time,
1308 *   Vector<double>& a_solution) const
1309 *   {
1310 *  
1311 *   DataOutBase::VtkFlags vtk_flags;
1312 *   vtk_flags.time = a_time;
1313 *   vtk_flags.cycle = a_time_idx;
1314 *  
1315 *   DataOut<dim> data_out;
1316 *   data_out.set_flags(vtk_flags);
1317 *  
1318 *   data_out.attach_dof_handler(dof_handler);
1319 *   data_out.add_data_vector(a_solution, "U");
1320 *  
1321 *   data_out.build_patches();
1322 *  
1323 *   const std::string filename = "solution-"
1324 *   + Utilities::int_to_string(a_time_idx, 3) +
1325 *   ".vtk";
1326 *   std::ofstream output(filename.c_str());
1327 *   data_out.write_vtk(output);
1328 *   }
1329 *  
1330 * @endcode
1331 *
1332 * We define the geometry here, this is called on each processor
1333 * and doesn't change in time. Once doing AMR, this won't need
1334 * to exist anymore.
1335 *
1336 * @code
1337 *   template <int dim>
1338 *   void HeatEquation<dim>::define()
1339 *   {
1340 *   const unsigned int initial_global_refinement = 6;
1341 *  
1343 *   triangulation.refine_global (initial_global_refinement);
1344 *  
1345 *   setup_system();
1346 *  
1347 *   tmp.reinit (dof_handler.n_dofs());
1348 *   forcing_terms.reinit (dof_handler.n_dofs());
1349 *   }
1350 *  
1351 * @endcode
1352 *
1353 * Here we advance the solution forward in time. This is done
1354 * the same way as in the loop in @ref step_26 "step-26"'s run function.
1355 *
1356 * @code
1357 *   template<int dim>
1358 *   void HeatEquation<dim>::step(Vector<double>& braid_data,
1359 *   double deltaT,
1360 *   double a_time,
1361 *   int a_time_idx)
1362 *   {
1363 *   a_time += deltaT;
1364 *   ++a_time_idx;
1365 *  
1366 *   mass_matrix.vmult(system_rhs, braid_data);
1367 *  
1368 *   laplace_matrix.vmult(tmp, braid_data);
1369 *  
1370 *   system_rhs.add(-(1 - theta) * deltaT, tmp);
1371 *  
1372 *   #if DO_MFG
1373 *   RightHandSideMFG<dim> rhs_function;
1374 *   #else
1375 *   RightHandSide<dim> rhs_function;
1376 *   #endif
1377 *   rhs_function.set_time(a_time);
1378 *   VectorTools::create_right_hand_side(dof_handler,
1379 *   QGauss<dim>(fe.degree+1),
1380 *   rhs_function,
1381 *   tmp);
1382 *  
1383 *   forcing_terms = tmp;
1384 *   forcing_terms *= deltaT * theta;
1385 *  
1386 *   rhs_function.set_time(a_time - deltaT);
1387 *   VectorTools::create_right_hand_side(dof_handler,
1388 *   QGauss<dim>(fe.degree+1),
1389 *   rhs_function,
1390 *   tmp);
1391 *  
1392 *   forcing_terms.add(deltaT * (1 - theta), tmp);
1393 *   system_rhs += forcing_terms;
1394 *  
1395 *   system_matrix.copy_from(mass_matrix);
1396 *   system_matrix.add(theta * deltaT, laplace_matrix);
1397 *  
1398 *   constraints.condense (system_matrix, system_rhs);
1399 *  
1400 *   {
1401 *   #if DO_MFG
1402 * @endcode
1403 *
1404 * If we are doing the method of manufactured solutions
1405 * then we set the boundary conditions to the exact solution.
1406 * Otherwise the boundary conditions are zero.
1407 *
1408 * @code
1409 *   ExactValuesMFG<dim> boundary_values_function;
1410 *   #else
1411 *   BoundaryValues<dim> boundary_values_function;
1412 *   #endif
1413 *   boundary_values_function.set_time(a_time);
1414 *  
1415 *   std::map<types::global_dof_index, double> boundary_values;
1416 *   VectorTools::interpolate_boundary_values(dof_handler,
1417 *   0,
1418 *   boundary_values_function,
1419 *   boundary_values);
1420 *  
1421 *   MatrixTools::apply_boundary_values(boundary_values,
1422 *   system_matrix,
1423 *   braid_data,
1424 *   system_rhs);
1425 *   }
1426 *  
1427 *   solve_time_step(braid_data);
1428 *   }
1429 *  
1430 *   template<int dim>
1431 *   int HeatEquation<dim>::size() const
1432 *   {
1433 *   return dof_handler.n_dofs();
1434 *   }
1435 *  
1436 * @endcode
1437 *
1438 * This function computes the error for the time step when doing
1439 * the method of manufactured solutions. First the exact values
1440 * is calculated, then the difference per cell is computed for
1441 * the various norms, and the error is computed. This is written
1442 * out to a pretty table.
1443 *
1444 * @code
1445 *   template<int dim> void
1446 *   HeatEquation<dim>::process_solution(double a_time,
1447 *   int a_index,
1448 *   const Vector<double>& a_vector)
1449 *   {
1450 * @endcode
1451 *
1452 * Compute the exact value for the manufactured solution case
1453 *
1454 * @code
1455 *   ExactValuesMFG<dim> exact_function;
1456 *   exact_function.set_time(a_time);
1457 *  
1458 *   Vector<double> difference_per_cell (triangulation.n_active_cells());
1459 *   VectorTools::integrate_difference(dof_handler,
1460 *   a_vector,
1461 *   exact_function,
1462 *   difference_per_cell,
1463 *   QGauss<dim>(fe.degree+1),
1464 *   VectorTools::L2_norm);
1465 *  
1466 *   const double L2_error = VectorTools::compute_global_error(triangulation,
1467 *   difference_per_cell,
1468 *   VectorTools::L2_norm);
1469 *  
1470 *   VectorTools::integrate_difference(dof_handler,
1471 *   a_vector,
1472 *   exact_function,
1473 *   difference_per_cell,
1474 *   QGauss<dim>(fe.degree+1),
1475 *   VectorTools::H1_seminorm);
1476 *  
1477 *   const double H1_error = VectorTools::compute_global_error(triangulation,
1478 *   difference_per_cell,
1479 *   VectorTools::H1_seminorm);
1480 *  
1481 *   const QTrapezoid<1> q_trapez;
1482 *   const QIterated<dim> q_iterated (q_trapez, 5);
1483 *   VectorTools::integrate_difference (dof_handler,
1484 *   a_vector,
1485 *   exact_function,
1486 *   difference_per_cell,
1487 *   q_iterated,
1488 *   VectorTools::Linfty_norm);
1489 *   const double Linfty_error = VectorTools::compute_global_error(triangulation,
1490 *   difference_per_cell,
1491 *   VectorTools::Linfty_norm);
1492 *  
1493 *   const unsigned int n_active_cells = triangulation.n_active_cells();
1494 *   const unsigned int n_dofs = dof_handler.n_dofs();
1495 *  
1496 *   pout() << "Cycle " << a_index << ':'
1497 *   << std::endl
1498 *   << " Number of active cells: "
1499 *   << n_active_cells
1500 *   << std::endl
1501 *   << " Number of degrees of freedom: "
1502 *   << n_dofs
1503 *   << std::endl;
1504 *  
1505 *   convergence_table.add_value("cycle", a_index);
1506 *   convergence_table.add_value("cells", n_active_cells);
1507 *   convergence_table.add_value("dofs", n_dofs);
1508 *   convergence_table.add_value("L2", L2_error);
1509 *   convergence_table.add_value("H1", H1_error);
1510 *   convergence_table.add_value("Linfty", Linfty_error);
1511 *  
1512 *   convergence_table.set_precision("L2", 3);
1513 *   convergence_table.set_precision("H1", 3);
1514 *   convergence_table.set_precision("Linfty", 3);
1515 *  
1516 *   convergence_table.set_scientific("L2", true);
1517 *   convergence_table.set_scientific("H1", true);
1518 *   convergence_table.set_scientific("Linfty", true);
1519 *  
1520 *   convergence_table.set_tex_caption("cells", "\\# cells");
1521 *   convergence_table.set_tex_caption("dofs", "\\# dofs");
1522 *   convergence_table.set_tex_caption("L2", "@fL^2@f-error");
1523 *   convergence_table.set_tex_caption("H1", "@fH^1@f-error");
1524 *   convergence_table.set_tex_caption("Linfty", "@fL^\\infty@f-error");
1525 *  
1526 *   convergence_table.set_tex_format("cells", "r");
1527 *   convergence_table.set_tex_format("dofs", "r");
1528 *  
1529 *   std::cout << std::endl;
1530 *   convergence_table.write_text(std::cout);
1531 *  
1532 *   std::ofstream error_table_file("tex-conv-table.tex");
1533 *   convergence_table.write_tex(error_table_file);
1534 *   }
1535 * @endcode
1536
1537
1538<a name="ann-src/Utilities.cc"></a>
1539<h1>Annotated version of src/Utilities.cc</h1>
1540 *
1541 *
1542 *
1543 *
1544 * @code
1545 *   /* -----------------------------------------------------------------------------
1546 *   *
1547 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1548 *   * Copyright (C) 2018 by Joshua Christopher
1549 *   *
1550 *   * This file is part of the deal.II code gallery.
1551 *   *
1552 *   * -----------------------------------------------------------------------------
1553 *   */
1554 *  
1555 *   #include "Utilities.hh"
1556 *  
1557 *   #include <string>
1558 *   #include <fstream>
1559 *  
1560 *   #include <mpi.h>
1561 *  
1562 *   int procID = 0;
1563 *  
1564 * @endcode
1565 *
1566 * The shared variables
1567 *
1568
1569 *
1570 *
1571 * @code
1572 *   static std::string s_pout_filename ;
1573 *   static std::string s_pout_basename ;
1574 *   static std::ofstream s_pout ;
1575 *  
1576 *   static bool s_pout_init = false ;
1577 *   static bool s_pout_open = false ;
1578 *  
1579 *   #ifdef USE_MPI
1580 * @endcode
1581 *
1582 * in parallel, compute the filename give the basename
1583 * [NOTE: dont call this before MPI is initialized.]
1584 *
1585 * @code
1586 *   static void setFileName()
1587 *   {
1588 *   static const size_t ProcnumSize = 1 + 10 + 1 ; //'.' + 10digits + '\0'
1589 *   char procnum[ProcnumSize] ;
1590 *   snprintf( procnum ,ProcnumSize ,".%d" ,procID);
1591 *   s_pout_filename = s_pout_basename + procnum ;
1592 *   }
1593 *  
1594 * @endcode
1595 *
1596 * in parallel, close the file if nec., open it and check for success
1597 *
1598 * @code
1599 *   static void openFile()
1600 *   {
1601 *   if ( s_pout_open )
1602 *   {
1603 *   s_pout.close();
1604 *   }
1605 *   s_pout.open( s_pout_filename.c_str() );
1606 * @endcode
1607 *
1608 * if open() fails, we have problems, but it's better
1609 * to try again later than to make believe it succeeded
1610 *
1611 * @code
1612 *   s_pout_open = (bool)s_pout ;
1613 *   }
1614 *  
1615 *   #else
1616 * @endcode
1617 *
1618 * in serial, filename is always cout
1619 *
1620 * @code
1621 *   static void setFileName()
1622 *   {
1623 *   s_pout_filename = "cout" ;
1624 *   }
1625 *  
1626 * @endcode
1627 *
1628 * in serial, this does absolutely nothing
1629 *
1630 * @code
1631 *   static void openFile()
1632 *   {
1633 *   }
1634 *   #endif
1635 *  
1636 *   std::ostream& pout()
1637 *   {
1638 *   #ifdef USE_MPI
1639 * @endcode
1640 *
1641 * the common case is _open == true, which just returns s_pout
1642 *
1643 * @code
1644 *   if ( ! s_pout_open )
1645 *   {
1646 * @endcode
1647 *
1648 * the uncommon cae: the file isn't opened, MPI may not be
1649 * initialized, and the basename may not have been set
1650 *
1651 * @code
1652 *   int flag_i, flag_f;
1653 *   MPI_Initialized(&flag_i);
1654 *   MPI_Finalized(&flag_f);
1655 * @endcode
1656 *
1657 * app hasn't set a basename yet, so set the default
1658 *
1659 * @code
1660 *   if ( ! s_pout_init )
1661 *   {
1662 *   s_pout_basename = "pout" ;
1663 *   s_pout_init = true ;
1664 *   }
1665 * @endcode
1666 *
1667 * if MPI not initialized, we cant open the file so return cout
1668 *
1669 * @code
1670 *   if ( ! flag_i || flag_f)
1671 *   {
1672 *   return std::cout; // MPI hasn't been started yet, or has ended....
1673 *   }
1674 * @endcode
1675 *
1676 * MPI is initialized, so file must not be, so open it
1677 *
1678 * @code
1679 *   setFileName() ;
1680 *   openFile() ;
1681 * @endcode
1682 *
1683 * finally, in case the open failed, return cout
1684 *
1685 * @code
1686 *   if ( ! s_pout_open )
1687 *   {
1688 *   return std::cout ;
1689 *   }
1690 *   }
1691 *   return s_pout ;
1692 *   #else
1693 *   return std::cout;
1694 *   #endif
1695 *   }
1696 * @endcode
1697
1698
1699<a name="ann-src/Utilities.hh"></a>
1700<h1>Annotated version of src/Utilities.hh</h1>
1701 *
1702 *
1703 *
1704 *
1705 * @code
1706 *   /* -----------------------------------------------------------------------------
1707 *   *
1708 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1709 *   * Copyright (C) 2018 by Joshua Christopher
1710 *   *
1711 *   * This file is part of the deal.II code gallery.
1712 *   *
1713 *   * -----------------------------------------------------------------------------
1714 *   */
1715 *  
1716 *   #ifndef _UTILITIES_H_
1717 *   #define _UTILITIES_H_
1718 *  
1719 *   #include <iostream>
1720 *  
1721 * @endcode
1722 *
1723 * This preprocessor macro is used on function arguments
1724 * that are not used in the function. It is used to
1725 * suppress compiler warnings.
1726 *
1727 * @code
1728 *   #define UNUSED(x) (void)(x)
1729 *  
1730 * @endcode
1731 *
1732 * Contains the current MPI processor ID.
1733 *
1734 * @code
1735 *   extern int procID;
1736 *  
1737 * @endcode
1738 *
1739 * Function to return the ostream to write out to. In MPI
1740 * mode it returns a stream to a file named pout.<#> where
1741 * <#> is the procID. This allows the user to write output
1742 * from each processor to a separate file. In serial mode
1743 * (no MPI), it returns the standard output.
1744 *
1745 * @code
1746 *   std::ostream& pout();
1747 *   #endif
1748 * @endcode
1749
1750
1751<a name="ann-src/parallel_in_time.cc"></a>
1752<h1>Annotated version of src/parallel_in_time.cc</h1>
1753 *
1754 *
1755 *
1756 *
1757 * @code
1758 *   /* -----------------------------------------------------------------------------
1759 *   *
1760 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1761 *   * Copyright (C) 2018 by Joshua Christopher
1762 *   *
1763 *   * This file is part of the deal.II code gallery.
1764 *   *
1765 *   * -----------------------------------------------------------------------------
1766 *   *
1767 *   * Author: Joshua Christopher, Colorado State University, 2018
1768 *   */
1769 *  
1770 *   #include "BraidFuncs.hh"
1771 *   #include "HeatEquation.hh"
1772 *   #include "Utilities.hh"
1773 *  
1774 *   #include <fstream>
1775 *   #include <iostream>
1776 *  
1777 *   int main(int argc, char *argv[])
1778 *   {
1779 *   try
1780 *   {
1781 *   using namespace dealii;
1782 *  
1783 *   /* Initialize MPI */
1784 *   MPI_Comm comm; //, comm_x, comm_t;
1785 *   int rank;
1786 *   MPI_Init(&argc, &argv);
1787 *   comm = MPI_COMM_WORLD;
1788 *   MPI_Comm_rank(comm, &rank);
1789 *   procID = rank;
1790 *  
1791 * @endcode
1792 *
1793 * Set up X-Braid
1794 *
1795 * @code
1796 *   /* Initialize Braid */
1797 *   braid_Core core;
1798 *   double tstart = 0.0;
1799 *   double tstop = 0.002;
1800 *   int ntime = 10;
1801 *   my_App *app = new(my_App);
1802 *  
1803 *   braid_Init(MPI_COMM_WORLD, comm, tstart, tstop, ntime, app,
1804 *   my_Step, my_Init, my_Clone, my_Free, my_Sum, my_SpatialNorm,
1805 *   my_Access, my_BufSize, my_BufPack, my_BufUnpack, &core);
1806 *  
1807 *   /* Define XBraid parameters
1808 *   * See -help message forf descriptions */
1809 *   int max_levels = 3;
1810 * @endcode
1811 *
1812 * int nrelax = 1;
1813 * int skip = 0;
1814 *
1815 * @code
1816 *   double tol = 1.e-7;
1817 * @endcode
1818 *
1819 * int cfactor = 2;
1820 *
1821 * @code
1822 *   int max_iter = 5;
1823 * @endcode
1824 *
1825 * int min_coarse = 10;
1826 * int fmg = 0;
1827 * int scoarsen = 0;
1828 * int res = 0;
1829 * int wrapper_tests = 0;
1830 *
1831 * @code
1832 *   int print_level = 1;
1833 *   int access_level = 1;
1834 *   int use_sequential= 0;
1835 *  
1836 *   braid_SetPrintLevel( core, print_level);
1837 *   braid_SetAccessLevel( core, access_level);
1838 *   braid_SetMaxLevels(core, max_levels);
1839 * @endcode
1840 *
1841 * braid_SetMinCoarse( core, min_coarse );
1842 * braid_SetSkip(core, skip);
1843 * braid_SetNRelax(core, -1, nrelax);
1844 *
1845 * @code
1846 *   braid_SetAbsTol(core, tol);
1847 * @endcode
1848 *
1849 * braid_SetCFactor(core, -1, cfactor);
1850 *
1851 * @code
1852 *   braid_SetMaxIter(core, max_iter);
1853 *   braid_SetSeqSoln(core, use_sequential);
1854 *  
1855 *   app->eq.define();
1856 *   app->final_step = ntime;
1857 *  
1858 *   braid_Drive(core);
1859 *  
1860 * @endcode
1861 *
1862 * Free the memory now that we are done
1863 *
1864 * @code
1865 *   braid_Destroy(core);
1866 *  
1867 *   delete app;
1868 *  
1869 * @endcode
1870 *
1871 * Clean up MPI
1872 * MPI_Comm_free(&comm);
1873 *
1874 * @code
1875 *   MPI_Finalize();
1876 *   }
1877 *   catch (std::exception &exc)
1878 *   {
1879 *   std::cerr << std::endl << std::endl
1880 *   << "----------------------------------------------------"
1881 *   << std::endl;
1882 *   std::cerr << "Exception on processing: " << std::endl << exc.what()
1883 *   << std::endl << "Aborting!" << std::endl
1884 *   << "----------------------------------------------------"
1885 *   << std::endl;
1886 *  
1887 *   return 1;
1888 *   }
1889 *   catch (...)
1890 *   {
1891 *   std::cerr << std::endl << std::endl
1892 *   << "----------------------------------------------------"
1893 *   << std::endl;
1894 *   std::cerr << "Unknown exception!" << std::endl << "Aborting!"
1895 *   << std::endl
1896 *   << "----------------------------------------------------"
1897 *   << std::endl;
1898 *   return 1;
1899 *   }
1900 *  
1901 *   return 0;
1902 *   }
1903 *  
1904 * @endcode
1905
1906
1907<a name="ann-test/test_braid.cc"></a>
1908<h1>Annotated version of test/test_braid.cc</h1>
1909 *
1910 *
1911 *
1912 *
1913 * @code
1914 *   /* -----------------------------------------------------------------------------
1915 *   *
1916 *   * SPDX-License-Identifier: LGPL-2.1-or-later
1917 *   * Copyright (C) 2018 by Joshua Christopher
1918 *   *
1919 *   * This file is part of the deal.II code gallery.
1920 *   *
1921 *   * -----------------------------------------------------------------------------
1922 *   */
1923 *  
1924 *   #include "BraidFuncs.hh"
1925 *  
1926 *   #include <braid.h>
1927 *   #include <braid_test.h>
1928 *  
1929 *   #include <iostream>
1930 *  
1931 *   int main(int argc, char** argv)
1932 *   {
1933 *   MPI_Comm comm;
1934 *   int rank;
1935 *   MPI_Init(&argc, &argv);
1936 *   comm = MPI_COMM_WORLD;
1937 *   MPI_Comm_rank(comm, &rank);
1938 *  
1939 *   my_App *app = new(my_App);
1940 *   app->eq.define();
1941 *  
1942 *   double time = 0.2;
1943 *  
1944 *   braid_Int init_access_result = braid_TestInitAccess(app,
1945 *   comm,
1946 *   stdout,
1947 *   time,
1948 *   my_Init,
1949 *   my_Access,
1950 *   my_Free);
1951 *   (void)init_access_result;
1952 *  
1953 *   braid_Int clone_result = braid_TestClone(app,
1954 *   comm,
1955 *   stdout,
1956 *   time,
1957 *   my_Init,
1958 *   my_Access,
1959 *   my_Free,
1960 *   my_Clone);
1961 *   (void)clone_result;
1962 *  
1963 *   braid_Int sum_result = braid_TestSum(app,
1964 *   comm,
1965 *   stdout,
1966 *   time,
1967 *   my_Init,
1968 *   my_Access,
1969 *   my_Free,
1970 *   my_Clone,
1971 *   my_Sum);
1972 *   (void)sum_result;
1973 *  
1974 *   braid_Int norm_result = braid_TestSpatialNorm(app,
1975 *   comm,
1976 *   stdout,
1977 *   time,
1978 *   my_Init,
1979 *   my_Free,
1980 *   my_Clone,
1981 *   my_Sum,
1982 *   my_SpatialNorm);
1983 *   (void)norm_result;
1984 *  
1985 *   braid_Int buf_result = braid_TestBuf(app,
1986 *   comm,
1987 *   stdout,
1988 *   time,
1989 *   my_Init,
1990 *   my_Free,
1991 *   my_Sum,
1992 *   my_SpatialNorm,
1993 *   my_BufSize,
1994 *   my_BufPack,
1995 *   my_BufUnpack);
1996 *   (void)buf_result;
1997 *  
1998 * @endcode
1999 *
2000 * /* Create spatial communicator for wrapper-tests */
2001 * braid_SplitCommworld(&comm, 1, &comm_x, &comm_t);
2002 *
2003 *
2004 * braid_TestAll(app, comm_x, stdout, 0.0, (tstop-tstart)/ntime,
2005 * 2*(tstop-tstart)/ntime, my_Init, my_Free, my_Clone,
2006 * my_Sum, my_SpatialNorm, my_BufSize, my_BufPack,
2007 * my_BufUnpack, my_Coarsen, my_Interp, my_Residual, my_Step);
2008 *
2009
2010 *
2011 *
2012 * @code
2013 *   /* Finalize MPI */
2015 *   }
2016 * @endcode
2017
2018
2019*/
void set_flags(const FlagType &flags)
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
unsigned int level
Definition grid_out.cc:4626
__global__ void set(Number *val, const Number s, const size_type N)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:442
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
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)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
const std::vector< bool > & used
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition l2.h:57
void create_mass_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, MatrixType &matrix, const Function< spacedim, typename MatrixType::value_type > *const a=nullptr, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >())
void create_laplace_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, MatrixType &matrix, const Function< spacedim, typename MatrixType::value_type > *const a=nullptr, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >())
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
const InputIterator OutputIterator out
Definition parallel.h:167
const InputIterator OutputIterator const Function & function
Definition parallel.h:168
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
*** braid_TestAll(app, comm_x, stdout, 0.0,(tstop-tstart)/ntime, *2 *(tstop-tstart)/ntime, my_Init, my_Free, my_Clone, *my_Sum, my_SpatialNorm, my_BufSize, my_BufPack, *my_BufUnpack, my_Coarsen, my_Interp, my_Residual, my_Step)
****code *  *  MPI_Finalize()
*braid_SplitCommworld & comm
void advance(std::tuple< I1, I2 > &t, const unsigned int n)