342 * braid_StepStatusGetLevel(status, &
level);
343 * braid_StepStatusGetTstartTstop(status, &tstart, &tstop);
344 * braid_StepStatusGetTIndex(status, &index);
346 * deltaT = tstop - tstart;
350 * HeatEquation<2>& heateq = app->eq;
352 * heateq.step(solution, deltaT, tstart, index);
360 * In
this function we initialize a vector at an arbitrary time.
361 * At
this point we don
't know anything about what the solution
362 * looks like, and we can really initialize to anything, so in
363 * this case use reinit to initialize the memory and set the
368 * my_Init(braid_App app,
370 * braid_Vector *u_ptr)
372 * my_Vector *u = new(my_Vector);
373 * int size = app->eq.size();
374 * u->data.reinit(size);
376 * app->eq.initialize(t, u->data);
385 * Here we need to copy the vector u into the vector v. We do this
386 * by allocating a new vector, then reinitializing the deal.ii
387 * vector to the correct size. The deal.ii reinitialization sets
388 * every value to zero, so next we need to iterate over the vector
389 * u and copy the values to the new vector v.
393 * my_Clone(braid_App app,
395 * braid_Vector *v_ptr)
398 * my_Vector *v = new(my_Vector);
399 * int size = u->data.size();
400 * v->data.reinit(size);
401 * for(size_t i=0, end=v->data.size(); i != end; ++i)
403 * v->data[i] = u->data[i];
412 * Here we need to free the memory used by vector u. This is
413 * pretty simple since the deal.ii vector is stored inside the
414 * XBraid vector, so we just delete the XBraid vector u and it
415 * puts the deal.ii vector out of scope and releases its memory.
419 * my_Free(braid_App app,
430 * This is to perform an axpy type operation. That is to say we
431 * do @f$y = \alpha x + \beta y@f$. Fortunately deal.ii already has
432 * this operation built in to its vector class, so we get the
433 * reference to the vector y and call the sadd method.
436 * int my_Sum(braid_App app,
443 * Vector<double>& vec = y->data;
444 * vec.sadd(beta, alpha, x->data);
451 * This calculates the spatial norm using the l2 norm. According
452 * to XBraid, this could be just about any spatial norm but we'll
453 * keep it simple and used deal.ii vector
's built in l2_norm method.
457 * my_SpatialNorm(braid_App app,
463 * dot = u->data.l2_norm();
471 * This function is called at various points depending on the access
472 * level specified when configuring the XBraid struct. This function
473 * is used to print out data during the run time, such as plots of the
474 * data. The status struct contains a ton of information about the
475 * simulation run. Here we get the current time and timestep number.
476 * The output_results function is called to plot the solution data.
477 * If the method of manufactured solutions is being used, then the
478 * error of this time step is computed and processed.
482 * my_Access(braid_App app,
484 * braid_AccessStatus astatus)
489 * braid_AccessStatusGetT(astatus, &t);
490 * braid_AccessStatusGetTIndex(astatus, &index);
492 * app->eq.output_results(index, t, u->data);
495 * if(index == app->final_step)
497 * app->eq.process_solution(t, index, u->data);
506 * This calculates the size of buffer needed to pack the solution
507 * data into a linear buffer for transfer to another processor via
508 * MPI. We query the size of the data from the HeatEquation class
509 * and return the buffer size.
513 * my_BufSize(braid_App app,
515 * braid_BufferStatus bstatus)
518 * int size = app->eq.size();
519 * *size_ptr = (size+1)*sizeof(double);
526 * This function packs a linear buffer with data so that the buffer
527 * may be sent to another processor via MPI. The buffer is cast to
528 * a type we can work with. The first element of the buffer is the
529 * size of the buffer. Then we iterate over soltuion vector u and
530 * fill the buffer with our solution data. Finally we tell XBraid
531 * how much data we wrote.
535 * my_BufPack(braid_App app,
538 * braid_BufferStatus bstatus)
542 * double *dbuffer = (double*)buffer;
543 * int size = u->data.size();
545 * for(int i=0; i != size; ++i)
547 * dbuffer[i+1] = (u->data)[i];
549 * braid_BufferStatusSetSize(bstatus, (size+1)*sizeof(double));
556 * This function unpacks a buffer that was recieved from a different
557 * processor via MPI. The size of the buffer is read from the first
558 * element, then we iterate over the size of the buffer and fill
559 * the values of solution vector u with the data in the buffer.
563 * my_BufUnpack(braid_App app,
565 * braid_Vector *u_ptr,
566 * braid_BufferStatus bstatus)
571 * my_Vector *u = NULL;
572 * double *dbuffer = (double*)buffer;
573 * int size = static_cast<int>(dbuffer[0]);
574 * u = new(my_Vector);
575 * u->data.reinit(size);
577 * for(int i = 0; i != size; ++i)
579 * (u->data)[i] = dbuffer[i+1];
588<a name="ann-src/BraidFuncs.hh"></a>
589<h1>Annotated version of src/BraidFuncs.hh</h1>
595 * #ifndef _BRAIDFUNCS_H_
596 * #define _BRAIDFUNCS_H_
599 * * \file BraidFuncs.cc
600 * * \brief Contains the implementation of the mandatory X-Braid functions
602 * * X-Braid mandates several functions in order to drive the solution.
603 * * This file contains the implementation of said mandatory functions.
604 * * See the X-Braid documentation for more information.
605 * * There are several functions that are optional in X-Braid that may
606 * * or may not be implemented in here.
611 * /*-------- Third Party --------*/
612 * #include <deal.II/numerics/vector_tools.h>
615 * #include <braid_test.h>
617 * /*-------- Project --------*/
618 * #include "HeatEquation.hh"
622 * This struct contains all data that changes with time. For now
623 * this is just the solution data. When doing AMR this should
624 * probably include the triangulization, the sparsity patter,
629 * * \brief Struct that contains the deal.ii vector.
631 * typedef struct _braid_Vector_struct
633 * ::Vector<double> data;
638 * This struct contains all the data that is unchanging with time.
642 * * \brief Struct that contains the HeatEquation and final
643 * * time step number.
645 * typedef struct _braid_App_struct
647 * HeatEquation<2> eq;
653 * * @brief my_Step - Takes a step in time, advancing the u vector
655 * * @param app - The braid app struct
656 * * @param ustop - The solution data at the end of this time step
657 * * @param fstop - RHS data (such as forcing function?)
658 * * @param u - The solution data at the beginning of this time step
659 * * @param status - Status structure that contains various info of this time
661 * * @return Success (0) or failure (1)
663 * int my_Step(braid_App app,
664 * braid_Vector ustop,
665 * braid_Vector fstop,
667 * braid_StepStatus status);
671 * * @brief my_Init - Initializes a solution data at the given time
672 * * For now, initializes the solution to zero no matter what time we are at
674 * * @param app - The braid app struct containing user data
675 * * @param t - Time at which the solution is initialized
676 * * @param u_ptr - The solution data that needs to be filled
678 * * @return Success (0) or failure (1)
681 * my_Init(braid_App app,
683 * braid_Vector *u_ptr);
687 * * @brief my_Clone - Clones a vector into a new vector
689 * * @param app - The braid app struct containing user data
690 * * @param u - The existing vector containing data
691 * * @param v_ptr - The empty vector that needs to be filled
693 * * @return Success (0) or failure (1)
696 * my_Clone(braid_App app,
698 * braid_Vector *v_ptr);
702 * * @brief my_Free - Deletes a vector
704 * * @param app - The braid app struct containing user data
705 * * @param u - The vector that needs to be deleted
707 * * @return Success (0) or failure (1)
710 * my_Free(braid_App app,
715 * * @brief my_Sum - Sums two vectors in an AXPY operation
716 * * The operation is y = alpha*x + beta*y
718 * * @param app - The braid app struct containing user data
719 * * @param alpha - The coefficient in front of x
720 * * @param x - A vector that is multiplied by alpha then added to y
721 * * @param beta - The coefficient of y
722 * * @param y - A vector that is multiplied by beta then summed with x
724 * * @return Success (0) or failure (1)
727 * my_Sum(braid_App app,
734 * * \brief Returns the spatial norm of the provided vector
736 * * Calculates and returns the spatial norm of the provided vector.
737 * * Interestingly enough, X-Braid does not specify a particular norm.
738 * * to keep things simple, we implement the Euclidean norm.
740 * * \param app - The braid app struct containing user data
741 * * \param u - The vector we need to take the norm of
742 * * \param norm_ptr - Pointer to the norm that was calculated, need to modify this
743 * * \return Success (0) or failure (1)
746 * my_SpatialNorm(braid_App app,
751 * * \brief Allows the user to output details
753 * * The Access function is called at various points to allow the user to output
754 * * information to the screen or to files.
755 * * The astatus parameter provides various information about the simulation,
756 * * see the XBraid documentation for details on what information you can get.
757 * * Example information is what the current timestep number and current time is.
758 * * If the access level (in parallel_in_time.cc) is set to 0, this function is
760 * * If the access level is set to 1, the function is called after the last
762 * * If the access level is set to 2, it is called every XBraid cycle.
764 * * \param app - The braid app struct containing user data
765 * * \param u - The vector containing the data at the status provided
766 * * \param astatus - The Braid status structure
767 * * \return Success (0) or failure (1)
770 * my_Access(braid_App app,
772 * braid_AccessStatus astatus);
775 * * \brief Calculates the size of a buffer for MPI data transfer
777 * * Calculates the size of the buffer that is needed to transfer
778 * * a solution vector to another processor.
779 * * The bstatus parameter provides various information on the
780 * * simulation, see the XBraid documentation for all possible
783 * * \param app - The braid app struct containing user data
784 * * \param size_ptr A pointer to the calculated size
785 * * \param bstatus The XBraid status structure
786 * * \return Success (0) or failure (1)
789 * my_BufSize(braid_App app,
791 * braid_BufferStatus bstatus);
794 * * \brief Linearizes a vector to be sent to another processor
796 * * Linearizes (packs) a data buffer with the contents of
797 * * some solution state u.
799 * * \param app - The braid app struct containing user data
800 * * \param u The vector that must be packed into buffer
801 * * \param buffer The buffer that must be filled with u
802 * * \param bstatus The XBraid status structure
803 * * \return Success (0) or failure (1)
806 * my_BufPack(braid_App app,
809 * braid_BufferStatus bstatus);
812 * * \brief Unpacks a vector that was sent from another processor
814 * * Unpacks a linear data buffer into the vector pointed to by
817 * * \param app - The braid app struct containing user data
818 * * \param buffer The buffer that must be unpacked
819 * * \param u_ptr The pointer to the vector that is filled
820 * * \param bstatus The XBraid status structure
821 * * \return Success (0) or failure (1)
824 * my_BufUnpack(braid_App app,
826 * braid_Vector *u_ptr,
827 * braid_BufferStatus bstatus);
829 * #endif // _BRAIDFUNCS_H_
833<a name="ann-src/HeatEquation.hh"></a>
834<h1>Annotated version of src/HeatEquation.hh</h1>
840 * #ifndef _HEATEQUATION_H_
841 * #define _HEATEQUATION_H_
843 * #include <deal.II/base/utilities.h>
844 * #include <deal.II/base/quadrature_lib.h>
845 * #include <deal.II/base/function.h>
846 * #include <deal.II/base/logstream.h>
847 * #include <deal.II/lac/vector.h>
848 * #include <deal.II/lac/full_matrix.h>
849 * #include <deal.II/lac/dynamic_sparsity_pattern.h>
850 * #include <deal.II/lac/sparse_matrix.h>
851 * #include <deal.II/lac/solver_cg.h>
852 * #include <deal.II/lac/precondition.h>
853 * #include <deal.II/lac/affine_constraints.h>
854 * #include <deal.II/grid/tria.h>
855 * #include <deal.II/grid/grid_generator.h>
856 * #include <deal.II/grid/grid_refinement.h>
857 * #include <deal.II/grid/grid_out.h>
858 * #include <deal.II/grid/tria_accessor.h>
859 * #include <deal.II/grid/tria_iterator.h>
860 * #include <deal.II/dofs/dof_handler.h>
861 * #include <deal.II/dofs/dof_accessor.h>
862 * #include <deal.II/dofs/dof_tools.h>
863 * #include <deal.II/fe/fe_q.h>
864 * #include <deal.II/fe/fe_values.h>
865 * #include <deal.II/numerics/data_out.h>
866 * #include <deal.II/numerics/vector_tools.h>
867 * #include <deal.II/numerics/error_estimator.h>
868 * #include <deal.II/numerics/solution_transfer.h>
869 * #include <deal.II/numerics/matrix_tools.h>
870 * #include <deal.II/base/convergence_table.h>
874 * using namespace dealii;
878 * The HeatEquation class is describes the finite element
879 * solver for the heat equation. It contains all the functions
880 * needed to define the problem domain and advance the solution
890 * void step(Vector<double>& braid_data,
895 * int size() const; /// Returns the size of the solution vector
897 * void output_results(int a_time_idx,
899 * Vector<double>& a_solution) const;
901 * void initialize(double a_time,
902 * Vector<double>& a_vector) const;
904 * void process_solution(double a_time,
906 * const Vector<double>& a_vector);
909 * void setup_system();
910 * void solve_time_step(Vector<double>& a_solution);
912 * Triangulation<dim> triangulation;
914 * DoFHandler<dim> dof_handler;
916 * AffineConstraints<double> constraints;
918 * SparsityPattern sparsity_pattern;
919 * SparseMatrix<double> mass_matrix;
920 * SparseMatrix<double> laplace_matrix;
921 * SparseMatrix<double> system_matrix;
923 * Vector<double> system_rhs;
925 * std::ofstream myfile;
927 * const double theta;
931 * These were originally in the run() function but because
932 * I am splitting the run() function up into define and step
933 * they need to become member data
936 * Vector<double> tmp;
937 * Vector<double> forcing_terms;
939 * ConvergenceTable convergence_table;
944 * The RightHandSide class describes the RHS of the governing
945 * equations. In this case, it is the forcing function.
949 * class RightHandSide : public Function<dim>
958 * virtual double value (const Point<dim> &p,
959 * const unsigned int component = 0) const;
962 * const double period;
967 * The BoundaryValues class describes the boundary conditions
968 * of the governing equations.
972 * class BoundaryValues : public Function<dim>
975 * virtual double value (const Point<dim> &p,
976 * const unsigned int component = 0) const;
981 * The RightHandSideMFG class describes the right hand side
982 * function when doing the method of manufactured solutions.
986 * class RightHandSideMFG : public Function<dim>
989 * virtual double value (const Point<dim> &p,
990 * const unsigned int component = 0) const;
995 * The InitialValuesMFG class describes the initial values
996 * when doing the method of manufactured solutions.
1000 * class InitialValuesMFG : public Function<dim>
1003 * virtual double value (const Point<dim> &p,
1004 * const unsigned int component = 0) const;
1009 * Provides the exact value for the manufactured solution. This
1010 * is used for the boundary conditions as well.
1013 * template <int dim>
1014 * class ExactValuesMFG : public Function<dim>
1018 * * \brief Computes the value at the given point and member data time
1020 * * Computes the exact value of the manufactured solution at point p and
1021 * * the member data time. See the class documentation and the design doc
1022 * * for details on what the exact solution is.
1024 * * \param p The point that the exact solution is computed at
1025 * * \param component The component of the exact solution (always 0 for now)
1026 * * \return double The exact value that was computed
1028 * virtual double value (const Point<dim> &p,
1029 * const unsigned int component = 0) const;
1032 * * \brief Computes the gradient of the exact solution at the given point
1034 * * Computes the gradient of the exact/manufactured solution value at
1035 * * point p and member data time. See the design doc for details on
1036 * * what the gradient of the exact solution is
1038 * * \param p The point that the gradient is calculated at
1039 * * \param component The component of the system of equations this gradient is for
1040 * * \return Tensor<1,dim> A rank 1 tensor that contains the gradient
1041 * * in each spatial dimension
1043 * virtual Tensor<1,dim> gradient (const Point<dim> &p,
1044 * const unsigned int component = 0) const;
1048 * #include "HeatEquationImplem.hh"
1050 * #endif // _HEATEQUATION_H_
1054<a name="ann-src/HeatEquationImplem.hh"></a>
1055<h1>Annotated version of src/HeatEquationImplem.hh</h1>
1061 * #include "Utilities.hh"
1063 * #include <iomanip>
1068 * Calculates the forcing function for the RightHandSide. See the
1069 * documentation for the math.
1072 * template <int dim>
1073 * double RightHandSide<dim>::value (const Point<dim> &p,
1074 * const unsigned int component) const
1077 * Assert (component == 0, ExcIndexRange(component, 0, 1));
1078 * Assert (dim == 2, ExcNotImplemented());
1080 * double time = this->get_time();
1082 * if ((p[0] > 0.5) && (p[1] > -0.5))
1084 * return std::exp(-0.5*(time-0.125)*(time-0.125)/(0.005));
1086 * else if ((p[0] > -0.5) && (p[1] > 0.5))
1088 * return std::exp(-0.5*(time-0.375)*(time-0.375)/(0.005));
1095 * return 0; // No forcing function
1100 * Calculates the forcing function for the method of manufactured
1101 * solutions. See the documentation for the math.
1104 * template <int dim>
1105 * double RightHandSideMFG<dim>::value (const Point<dim> &p,
1106 * const unsigned int component) const
1109 * Assert (component == 0, ExcIndexRange(component, 0, 1));
1110 * Assert (dim == 2, ExcNotImplemented());
1112 * double time = this->get_time();
1114 * double pi = numbers::PI;
1115 * return 4*pi*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
1120 * Calculates the boundary conditions, essentially zero everywhere.
1123 * template <int dim>
1124 * double BoundaryValues<dim>::value (const Point<dim> &p,
1125 * const unsigned int component) const
1129 * Assert (component == 0, ExcIndexRange(component, 0, 1));
1135 * Calculates the exact solution (and thus also boundary conditions)
1136 * for the method of manufactured solutions.
1139 * template <int dim>
1140 * double ExactValuesMFG<dim>::value (const Point<dim> &p,
1141 * const unsigned int component) const
1144 * Assert (component == 0, ExcIndexRange(component, 0, 1));
1146 * double time = this->get_time();
1147 * const double pi = numbers::PI;
1149 * return std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
1154 * Calculates the gradient of the exact solution for the method of manufactured
1155 * solutions. See the documentation for the math.
1158 * template <int dim>
1159 * Tensor<1,dim> ExactValuesMFG<dim>::gradient (const Point<dim> &p,
1160 * const unsigned int) const
1162 * Assert (dim == 2, ExcNotImplemented());
1164 * Tensor<1,dim> return_value;
1165 * const double pi = numbers::PI;
1166 * double time = this->get_time();
1167 * return_value[0] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[1])*std::sin(2*pi*p[0]);
1168 * return_value[1] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::sin(2*pi*p[1]);
1169 * return return_value;
1174 * Calculates the initial values for the method of manufactured solutions.
1175 * See the documentation for the math.
1178 * template <int dim>
1179 * double InitialValuesMFG<dim>::value (const Point<dim> &p,
1180 * const unsigned int component) const
1183 * Assert (component == 0, ExcIndexRange(component, 0, 1));
1184 * const double pi = numbers::PI;
1186 * return std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
1189 * template <int dim>
1190 * HeatEquation<dim>::HeatEquation ()
1193 * dof_handler(triangulation),
1198 * template <int dim>
1199 * void HeatEquation<dim>::initialize(double a_time,
1200 * Vector<double>& a_vector) const
1205 * We only initialize values in the manufactured solution case
1208 * InitialValuesMFG<dim> iv_function;
1209 * iv_function.set_time(a_time);
1210 * VectorTools::project (dof_handler, constraints,
1211 * QGauss<dim>(fe.degree+1), iv_function,
1219 * If not the MFG solution case, a_vector is already zero'd so
do nothing
1224 *
template <
int dim>
1225 *
void HeatEquation<dim>::setup_system()
1227 * dof_handler.distribute_dofs(fe);
1229 * constraints.clear ();
1232 * constraints.close();
1239 * sparsity_pattern.copy_from(dsp);
1242 * laplace_matrix.reinit(sparsity_pattern);
1243 * system_matrix.reinit(sparsity_pattern);
1252 * system_rhs.reinit(dof_handler.n_dofs());
1256 *
template <
int dim>
1257 *
void HeatEquation<dim>::solve_time_step(
Vector<double>& a_solution)
1259 *
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
1263 * preconditioner.
initialize(system_matrix, 1.0);
1265 * cg.solve(system_matrix, a_solution, system_rhs,
1268 * constraints.distribute(a_solution);
1273 *
template <
int dim>
1274 *
void HeatEquation<dim>::output_results(
int a_time_idx,
1280 * vtk_flags.
time = a_time;
1281 * vtk_flags.
cycle = a_time_idx;
1291 *
const std::string filename =
"solution-"
1294 * std::ofstream output(filename.c_str());
1300 * We define the geometry here,
this is called on each processor
1301 * and doesn
't change in time. Once doing AMR, this won't need
1305 *
template <
int dim>
1306 *
void HeatEquation<dim>::define()
1308 *
const unsigned int initial_global_refinement = 6;
1315 * tmp.reinit (dof_handler.n_dofs());
1316 * forcing_terms.reinit (dof_handler.n_dofs());
1321 * Here we
advance the solution forward in time. This is done
1322 * the same way as in the
loop in @ref step_26
"step-26"'s run function.
1326 * void HeatEquation<dim>::step(Vector<double>& braid_data,
1334 * mass_matrix.vmult(system_rhs, braid_data);
1336 * laplace_matrix.vmult(tmp, braid_data);
1338 * system_rhs.add(-(1 - theta) * deltaT, tmp);
1341 * RightHandSideMFG<dim> rhs_function;
1343 * RightHandSide<dim> rhs_function;
1345 * rhs_function.set_time(a_time);
1346 * VectorTools::create_right_hand_side(dof_handler,
1347 * QGauss<dim>(fe.degree+1),
1351 * forcing_terms = tmp;
1352 * forcing_terms *= deltaT * theta;
1354 * rhs_function.set_time(a_time - deltaT);
1355 * VectorTools::create_right_hand_side(dof_handler,
1356 * QGauss<dim>(fe.degree+1),
1360 * forcing_terms.add(deltaT * (1 - theta), tmp);
1361 * system_rhs += forcing_terms;
1363 * system_matrix.copy_from(mass_matrix);
1364 * system_matrix.add(theta * deltaT, laplace_matrix);
1366 * constraints.condense (system_matrix, system_rhs);
1372 * If we are doing the method of manufactured solutions
1373 * then we set the boundary conditions to the exact solution.
1374 * Otherwise the boundary conditions are zero.
1377 * ExactValuesMFG<dim> boundary_values_function;
1379 * BoundaryValues<dim> boundary_values_function;
1381 * boundary_values_function.set_time(a_time);
1383 * std::map<types::global_dof_index, double> boundary_values;
1384 * VectorTools::interpolate_boundary_values(dof_handler,
1386 * boundary_values_function,
1389 * MatrixTools::apply_boundary_values(boundary_values,
1395 * solve_time_step(braid_data);
1399 * int HeatEquation<dim>::size() const
1401 * return dof_handler.n_dofs();
1406 * This function computes the error for the time step when doing
1407 * the method of manufactured solutions. First the exact values
1408 * is calculated, then the difference per cell is computed for
1409 * the various norms, and the error is computed. This is written
1410 * out to a pretty table.
1413 * template<int dim> void
1414 * HeatEquation<dim>::process_solution(double a_time,
1416 * const Vector<double>& a_vector)
1420 * Compute the exact value for the manufactured solution case
1423 * ExactValuesMFG<dim> exact_function;
1424 * exact_function.set_time(a_time);
1426 * Vector<double> difference_per_cell (triangulation.n_active_cells());
1427 * VectorTools::integrate_difference(dof_handler,
1430 * difference_per_cell,
1431 * QGauss<dim>(fe.degree+1),
1432 * VectorTools::L2_norm);
1434 * const double L2_error = VectorTools::compute_global_error(triangulation,
1435 * difference_per_cell,
1436 * VectorTools::L2_norm);
1438 * VectorTools::integrate_difference(dof_handler,
1441 * difference_per_cell,
1442 * QGauss<dim>(fe.degree+1),
1443 * VectorTools::H1_seminorm);
1445 * const double H1_error = VectorTools::compute_global_error(triangulation,
1446 * difference_per_cell,
1447 * VectorTools::H1_seminorm);
1449 * const QTrapez<1> q_trapez;
1450 * const QIterated<dim> q_iterated (q_trapez, 5);
1451 * VectorTools::integrate_difference (dof_handler,
1454 * difference_per_cell,
1456 * VectorTools::Linfty_norm);
1457 * const double Linfty_error = VectorTools::compute_global_error(triangulation,
1458 * difference_per_cell,
1459 * VectorTools::Linfty_norm);
1461 * const unsigned int n_active_cells = triangulation.n_active_cells();
1462 * const unsigned int n_dofs = dof_handler.n_dofs();
1464 * pout() << "Cycle " << a_index << ':
'
1466 * << " Number of active cells: "
1469 * << " Number of degrees of freedom: "
1473 * convergence_table.add_value("cycle", a_index);
1474 * convergence_table.add_value("cells", n_active_cells);
1475 * convergence_table.add_value("dofs", n_dofs);
1476 * convergence_table.add_value("L2", L2_error);
1477 * convergence_table.add_value("H1", H1_error);
1478 * convergence_table.add_value("Linfty", Linfty_error);
1480 * convergence_table.set_precision("L2", 3);
1481 * convergence_table.set_precision("H1", 3);
1482 * convergence_table.set_precision("Linfty", 3);
1484 * convergence_table.set_scientific("L2", true);
1485 * convergence_table.set_scientific("H1", true);
1486 * convergence_table.set_scientific("Linfty", true);
1488 * convergence_table.set_tex_caption("cells", "\\# cells");
1489 * convergence_table.set_tex_caption("dofs", "\\# dofs");
1490 * convergence_table.set_tex_caption("L2", "@fL^2@f-error");
1491 * convergence_table.set_tex_caption("H1", "@fH^1@f-error");
1492 * convergence_table.set_tex_caption("Linfty", "@fL^\\infty@f-error");
1494 * convergence_table.set_tex_format("cells", "r");
1495 * convergence_table.set_tex_format("dofs", "r");
1497 * std::cout << std::endl;
1498 * convergence_table.write_text(std::cout);
1500 * std::ofstream error_table_file("tex-conv-table.tex");
1501 * convergence_table.write_tex(error_table_file);
1506<a name="ann-src/Utilities.cc"></a>
1507<h1>Annotated version of src/Utilities.cc</h1>
1513 * #include "Utilities.hh"
1516 * #include <fstream>
1524 * The shared variables
1530 * static std::string s_pout_filename ;
1531 * static std::string s_pout_basename ;
1532 * static std::ofstream s_pout ;
1534 * static bool s_pout_init = false ;
1535 * static bool s_pout_open = false ;
1540 * in parallel, compute the filename give the basename
1541 * [NOTE: dont call this before MPI is initialized.]
1544 * static void setFileName()
1546 * static const size_t ProcnumSize = 1 + 10 + 1 ; //'.
' + 10digits + '\0
'
1547 * char procnum[ProcnumSize] ;
1548 * snprintf( procnum ,ProcnumSize ,".%d" ,procID);
1549 * s_pout_filename = s_pout_basename + procnum ;
1554 * in parallel, close the file if nec., open it and check for success
1557 * static void openFile()
1559 * if ( s_pout_open )
1563 * s_pout.open( s_pout_filename.c_str() );
1566 * if open() fails, we have problems, but it's better
1567 * to
try again later than to make believe it succeeded
1570 * s_pout_open = (
bool)s_pout ;
1579 *
static void setFileName()
1581 * s_pout_filename =
"cout" ;
1589 *
static void openFile()
1594 * std::ostream& pout()
1599 * the common
case is _open ==
true, which just returns s_pout
1602 *
if ( ! s_pout_open )
1606 * the uncommon cae: the file isn
't opened, MPI may not be
1607 * initialized, and the basename may not have been set
1610 * int flag_i, flag_f;
1611 * MPI_Initialized(&flag_i);
1612 * MPI_Finalized(&flag_f);
1615 * app hasn't
set a basename yet, so
set the
default
1618 *
if ( ! s_pout_init )
1620 * s_pout_basename =
"pout" ;
1621 * s_pout_init = true ;
1625 *
if MPI not initialized, we cant open the file so
return cout
1628 *
if ( ! flag_i || flag_f)
1634 * MPI is initialized, so file must not be, so open it
1641 *
finally, in
case the open failed,
return cout
1644 *
if ( ! s_pout_open )
1646 *
return std::cout ;
1657<a name=
"ann-src/Utilities.hh"></a>
1658<h1>Annotated version of src/
Utilities.hh</h1>
1664 * #ifndef _UTILITIES_H_
1665 * #define _UTILITIES_H_
1667 * #include <iostream>
1671 * This preprocessor macro is used on function arguments
1672 * that are not used in the function. It is used to
1673 * suppress compiler warnings.
1676 * #define UNUSED(x) (void)(x)
1680 * Contains the current MPI processor ID.
1683 *
extern int procID;
1687 *
Function to
return the ostream to write out to. In MPI
1688 * mode it returns a stream to a file named pout.<#> where
1689 * <#> is the procID. This allows the user to write output
1690 * from each processor to a separate file. In
serial mode
1691 * (no MPI), it returns the standard output.
1694 * std::ostream& pout();
1699<a name=
"ann-src/parallel_in_time.cc"></a>
1700<h1>Annotated version of src/parallel_in_time.cc</h1>
1725 * #include
"BraidFuncs.hh"
1726 * #include
"HeatEquation.hh"
1727 * #include
"Utilities.hh"
1729 * #include <fstream>
1730 * #include <iostream>
1732 *
int main(
int argc,
char *argv[])
1736 *
using namespace dealii;
1741 * MPI_Init(&argc, &argv);
1742 *
comm = MPI_COMM_WORLD;
1743 * MPI_Comm_rank(
comm, &rank);
1753 *
double tstart = 0.0;
1754 *
double tstop = 0.002;
1756 * my_App *app =
new(my_App);
1758 * braid_Init(MPI_COMM_WORLD,
comm, tstart, tstop, ntime, app,
1759 * my_Step, my_Init, my_Clone, my_Free, my_Sum, my_SpatialNorm,
1760 * my_Access, my_BufSize, my_BufPack, my_BufUnpack, &core);
1764 *
int max_levels = 3;
1771 *
double tol = 1.e-7;
1780 *
int min_coarse = 10;
1784 *
int wrapper_tests = 0;
1787 *
int print_level = 1;
1788 *
int access_level = 1;
1789 *
int use_sequential= 0;
1791 * braid_SetPrintLevel( core, print_level);
1792 * braid_SetAccessLevel( core, access_level);
1793 * braid_SetMaxLevels(core, max_levels);
1796 * braid_SetMinCoarse( core, min_coarse );
1797 * braid_SetSkip(core, skip);
1798 * braid_SetNRelax(core, -1, nrelax);
1801 * braid_SetAbsTol(core, tol);
1804 * braid_SetCFactor(core, -1, cfactor);
1807 * braid_SetMaxIter(core, max_iter);
1808 * braid_SetSeqSoln(core, use_sequential);
1811 * app->final_step = ntime;
1813 * braid_Drive(core);
1817 * Free the memory now that we are done
1820 * braid_Destroy(core);
1827 * MPI_Comm_free(&
comm);
1832 *
catch (std::exception &exc)
1834 * std::cerr << std::endl << std::endl
1835 * <<
"----------------------------------------------------"
1837 * std::cerr <<
"Exception on processing: " << std::endl << exc.what()
1838 * << std::endl <<
"Aborting!" << std::endl
1839 * <<
"----------------------------------------------------"
1846 * std::cerr << std::endl << std::endl
1847 * <<
"----------------------------------------------------"
1849 * std::cerr <<
"Unknown exception!" << std::endl <<
"Aborting!"
1851 * <<
"----------------------------------------------------"
1862<a name=
"ann-test/test_braid.cc"></a>
1863<h1>Annotated version of test/test_braid.cc</h1>
1869 * #include
"BraidFuncs.hh"
1871 * #include <braid.h>
1872 * #include <braid_test.h>
1874 * #include <iostream>
1876 *
int main(
int argc,
char** argv)
1880 * MPI_Init(&argc, &argv);
1881 *
comm = MPI_COMM_WORLD;
1882 * MPI_Comm_rank(
comm, &rank);
1884 * my_App *app =
new(my_App);
1887 *
double time = 0.2;
1889 * braid_Int init_access_result = braid_TestInitAccess(app,
1896 * (void)init_access_result;
1898 * braid_Int clone_result = braid_TestClone(app,
1906 * (void)clone_result;
1908 * braid_Int sum_result = braid_TestSum(app,
1919 * braid_Int norm_result = braid_TestSpatialNorm(app,
1928 * (void)norm_result;
1930 * braid_Int buf_result = braid_TestBuf(app,
1946 * braid_SplitCommworld(&
comm, 1, &comm_x, &comm_t);
1950 * 2*(tstop-tstart)/ntime, my_Init, my_Free, my_Clone,
1951 * my_Sum, my_SpatialNorm, my_BufSize, my_BufPack,
1952 * my_BufUnpack, my_Coarsen, my_Interp, my_Residual, my_Step);
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
__global__ void set(Number *val, const Number s, const size_type N)
void set_flags(const FlagType &flags)
void write_vtk(std::ostream &out) const
void loop(ITERATOR begin, typename identity< ITERATOR >::type 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())
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &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.)
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)
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)
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)