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