--- /dev/null
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+set (CMAKE_CXX_STANDARD 14)
+
+# Set the name of the project and target:
+SET(TARGET "parallel_in_time")
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -pedantic")
+
+
+set(PROJECT_PATH "${CMAKE_CURRENT_SOURCE_DIR}")
+set(SRC_PATH "${PROJECT_PATH}/src")
+set(TEST_SRC_PATH "${PROJECT_PATH}/test")
+# set(BIN_PATH "${PROJECT_BINARY_DIR}/bin")
+set(BIN_NAME "${TARGET}")
+set(TEST_PATH "${PROJECT_BINARY_DIR}/test")
+set(TEST_NAME "test_${TARGET}")
+
+SET(TARGET_SRC
+ ${SRC_PATH}/${TARGET}.cc
+ ${SRC_PATH}/BraidFuncs.cc
+ ${SRC_PATH}/Utilities.cc
+ )
+
+SET(TEST_SRC
+ ${TEST_SRC_PATH}/test_braid.cc
+ )
+
+OPTION(DO_MFG "Do the manufactured solution run" OFF)
+
+# Tell CMake to find MPI
+find_package(MPI REQUIRED)
+include_directories(${MPI_INCLUDE_PATH})
+add_definitions(-DUSE_MPI)
+
+if(DO_MFG)
+ # If we are doing the manufactured solution, tell
+ # the C++ processor so it can enable it
+ add_definitions(-DDO_MFG)
+endif(DO_MFG)
+
+FIND_PACKAGE(deal.II 9.0.0 QUIET
+ HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+ )
+IF(NOT ${deal.II_FOUND})
+ MESSAGE(FATAL_ERROR "\n"
+ "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+ "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+ "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+ )
+ENDIF()
+# include_directories(${DEAL_II_INCLUDE_DIRS})
+
+
+
+# Find braid details
+if (NOT "$ENV{BRAID_DIR}" STREQUAL "")
+ set(BRAID_DIR "$ENV{BRAID_DIR}" CACHE INTERNAL "Copied BRAID_DIR from environment variable")
+endif()
+# Include the braid paths and libraries
+include_directories(${BRAID_DIR})
+
+
+##############################
+# Finally start building stuff
+
+# First, build the "library" that consists of all the source files
+# set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${LIB_PATH})
+# set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${LIB_PATH})
+# add_library(${LIB_NAME} ${SRC})
+
+
+# Next build the main function
+# set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${BIN_PATH})
+# add_executable(${BIN_NAME} ${MAIN_SRC})
+# target_link_libraries(${MAIN_NAME} ${LIB_NAME}) # link my library
+# target_link_libraries(${MAIN_NAME} ${BRAID_DIR}libbraid.a) # link braid
+# DEAL_II_SETUP_TARGET(${MAIN_NAME}) # Link dealii
+# if(USE_MPI)
+# set_target_properties(${MAIN_NAME} PROPERTIES
+# LINK_FLAGS "${MPI_LINK_FLAGS}") # Link MPI
+# set_target_properties(${MAIN_NAME} PROPERTIES
+# COMPILE_FLAGS "${MPI_COMPILE_FLAGS}") # Use MPI compile flags
+# endif(USE_MPI)
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
+target_link_libraries(${TARGET} ${BRAID_DIR}libbraid.a) # link braid
+set_target_properties(${TARGET} PROPERTIES
+ LINK_FLAGS "${MPI_LINK_FLAGS}") # Link MPI
+set_target_properties(${MAIN_NAME} PROPERTIES
+ COMPILE_FLAGS "${MPI_COMPILE_FLAGS}") # Use MPI compile flags
+
+# Lastly build the tests
+enable_testing()
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${TEST_PATH})
+add_subdirectory(test)
+
+
+
--- /dev/null
+# Parallel in Time deal.ii
+
+## Overview
+
+Over the last few years, the clock speed per processor core has stagnated.
+This stagnation has lead to the design of larger core counts in high performance computing machines.
+As a result of these developments, increased concurrency in numerical algorithms must be developed in order to take advantage of this architecture style.
+Perhaps the largest bottleneck in concurrency for time-dependent simulations is traditional time integration methods.
+Traditional time integration methods solve for the time domain sequentially, and as the spatial grid is refined a proportionally larger number of time steps must be taken to maintain accuracy and stability constraints.
+While solving the time domain sequentially with a traditional time integration method is an optimal algorithm of order $\mathcal{O}(n)$, the $n$ time steps are not solved concurrently.
+
+The goal of this project is to make use of the XBraid library from Lawrence Livermore National Laboratory to solve the time domain in parallel using multigrid-reduction-in-time techniques.
+The XBraid library is implemented in C and aims to be a non-intrusive method to implement parallel time marching methods into existing codes.
+
+## Implementation
+
+### XBraid introduction
+
+In order to use the XBraid library, several data structures and functions must be implemented and provided to the XBraid solver struct.
+The two required data structures are the <kbd>app</kbd> and <kbd>vector</kbd> structures.
+In general, the <kbd>app</kbd> struct contains the time independent data and the <kbd>vector</kbd> struct contains the time dependent data.
+For this initial example, the time independent data includes the mesh which is fixed for all time steps, and the time dependent data is the solution state vector.
+The functions tell XBraid how to perform operations on the data type used by your solver, in this case deal.ii uses the Vector<double> data type.
+These operations include how to initialize the data at a given time, how to sum the data, and how to pack and unpack linear buffers for transmission to other processors via MPI.
+The XBraid documentation should be read for a full list of functions that must be implemented and the details of what the function should do.
+The typical format is the function is called with arguments of the <kbd>app</kbd> struct, one or more <kbd>vector</kbd> structs, and a <kbd>status</kbd> struct that contains information on the current status of the XBraid simulation (the current multigrid iteration, the level the function is being called from, the time and timestep number, etc.).
+
+Perhaps the most important function is the <kbd>step</kbd> function.
+This function tells XBraid how to advance the solution forward in time from the initial to the final times given in the <kbd>status</kbd> struct.
+This method uses a traditional time integration method such as the fourth order explicit Runge Kutta method.
+
+### deal.ii details
+
+The solver used in this example is based off the heat equation solver from the [Step-26 Tutorial](https://dealii.org/developer/doxygen/deal.II/step_26.html .
+The <kbd>HeatEquation</kbd> class becomes member data to XBraid's <kbd>app</kbd> struct, and XBraid's <kbd>vector</kbd> struct becomes a wrapper for deal.ii's <kbd>Vector<double></kbd> data type.
+The <kbd>HeatEquation</kbd> class cannot simply be used as is though as it contains both time dependent and time independent member data.
+In order to simplify the problem the adaptive mesh refinement is removed.
+Theoretically XBraid is capable of working with adaptive mesh refinement and in fact contains support for time refinement (which is also not used for simplicity).
+All adaptive mesh refinement functionality is removed from the sovler.
+The time-dependent solution state vectors are also removed from the <kbd>HeatEquation</kbd> member data.
+That time-dependent data will be provided at each timestep by XBraid via the <kbd>vector</kbd> struct.
+
+### Math details
+
+In the default mode, this code solves the heat equation with,
+
+@f{align}
+\frac{\partial u}{\partial t} - \Delta u = f(\boldsymbol{x},t), \qquad \forall\boldsymbol{x}\in\Omega,t\in\left( 0,T \right),
+@f}
+with initial conditions,
+@f{align}
+u(\boldsymbol{x},0) = u_0(\boldsymbol{x}) = 0, \qquad \forall \boldsymbol{x}\in\Omega,
+@f}
+and dirichlet boundary conditions,
+@f{align}
+u(\boldsymbol{x},t) = g(\boldsymbol{x},t) = 0, \qquad \forall \boldsymbol{x}\in\partial\Omega,t\in\left( 0,T \right)
+@f}
+and forcing function,
+@f{align}
+ f(\mathbf x, t) =
+ \left\{
+ \begin{array}{ll}
+ \chi_1(t) & \text{if \(x>0.5\) and \(y>-0.5\)} \\
+ \chi_2(t) & \text{if \(x>-0.5\) and \(y>0.5\)}
+ \end{array}
+ \right.
+@f}
+with,
+@f{align}
+ \chi_1(t) = \exp\left(-0.5\frac{(t-0.125)^2}{0.005}\right)
+@f}
+and,
+@f{align}
+ \chi_2(t) = \exp\left(-0.5\frac{(t-0.375)^2}{0.005}\right)
+@f}
+for some time $t$.
+The forcing function is a Gaussian pulse in time that is centered around 0.125 time units for $\chi_1$ and 0.375 time units for $\chi_2$.
+A gaussian function was chosen because it is a continuous function in time and so the solution state can be compared bitwise with the serial-in-time solution.
+
+The method of manufactured solutions is used to test the correctness of the implementation.
+In the method of manufactured solutions, we create a solution $u_h$ to the heat equation, then compute the boundary conditions, initial conditions, and forcing functions required to generate that solution.
+The created solution used is,
+@f{align}
+u_h = \exp\left(-4\pi^2t\right)\cos(2 \pi x)\cos(2 \pi y), \qquad \forall \boldsymbol{x} \in \Omega \cup \partial\Omega
+@f}
+with derivatives,
+@f{align}
+\frac{\partial u}{\partial t} &= -4 \pi^2 e^{-4\pi^2t}\cos(2 \pi x)\cos(2 \pi y), \\
+-\Delta u &= 8 \pi^2 \exp\left(-4\pi^2t\right)\cos(2 \pi x)\cos(2 \pi y) \\
+\frac{\partial u}{\partial x} &= -2 \pi \exp\left(-4\pi^2t\right)\sin(2\pi x)\cos(2\pi y) \\
+\frac{\partial u}{\partial x} &= -2 \pi \exp\left(-4\pi^2t\right)\cos(2\pi x)\sin(2\pi y)
+@f}
+and therefore we specify the forcing term, initial conditions, and boundary conditions of the governing equations as,
+@f{align}
+f(\boldsymbol{x},t) &= 4 \pi^2 \exp\left(-4\pi^2t\right)\cos(2 \pi x)\cos(2 \pi y), &\forall\boldsymbol{x}\in\Omega,t\in\left( 0,T \right), \\
+u_0(\boldsymbol{x}) &= \cos(2 \pi x)\cos(2 \pi y), &\forall \boldsymbol{x}\in\Omega, \\
+g(\boldsymbol{x},t) &= \exp\left(-4\pi^2t\right)\cos(2 \pi x)\cos(2 \pi y), &\forall \boldsymbol{x} \in \partial\Omega.
+@f}
+
+The manufactured solution is run on progressively more refined grids and the solution generated by the finite element method is compared to the exact solution $u_h$.
+The convergence rate of the error is calculated with
+@f{align}
+ \Delta \epsilon_n = \frac{\ln{\epsilon_{n-1}/\epsilon_{n}}}{\ln{r_n}}
+@f}
+where $\Delta \epsilon_n$ is the convergence rate of error $\epsilon$ between a mesh $n$ and
+coarser mesh $n-1$ that have a refinement ratio of $r_n$.
+Shown in Table 1 is the convergence rate as the mesh is refined.
+The $\Delta t$ is reduced by a factor of 2 for every global refinement of the mesh.
+
+@f{align}
+\begin{tabular}{|c|c|c|c|c|c|} \hline
+ cycle & \# cells & \# dofs & L^2-error & H^1-error & L^\infty-error \\ \hline
+ 125 & 48 & 65 & 6.036e-03 & 6.970e-02 & 7.557e-03\\ \hline
+ 250 & 192 & 225 & 1.735e-03 & 3.414e-02 & 2.721e-03 \\ \hline
+ 500 & 768 & 833 & 4.513e-04 & 1.690e-02 & 7.410e-04 \\ \hline
+ 1000 & 3072 & 3201 & 1.140e-04 & 8.426e-03 & 1.877e-04 \\ \hline
+ 2000 & 12288 & 12545 & 2.859e-05 & 4.209e-03 & 4.715e-05 \\ \hline
+\end{tabular}
+@f}
+
+@f{align}
+\begin{tabular}{|c|c|c|c|c|c|} \hline
+ cycle & \# cells & \# dofs & Slope L^2 & Slope H^1 & Slope \textrm{L}^\infty \\ \hline
+ 125 & 48 & 65 & --- & --- & --- \\ \hline
+ 250 & 192 & 225 & 1.798 & 1.030 & 1.474 \\ \hline
+ 500 & 768 & 833 & 1.943 & 1.014 & 1.877 \\ \hline
+ 1000 & 3072 & 3201 & 1.985 & 1.004 & 1.981 \\ \hline
+ 2000 & 12288 & 12545 & 1.995 & 1.001 & 1.993 \\ \hline
+\end{tabular}
+@f}
+
+As can be seen, the slope converges at a second order rate for the $\textrm{L}^2$ and $\textrm{L}^\infty$ norms and at a first order rate for the $\textrm{H}^1$ norms.
+This is expected behavior as second order finite elements are used.
+
+### Code Organization
+
+#### The src directory
+
+The entry point of the code is in <kbd>parallel_in_time.cc</kbd> and sets up XBraid for a simulation. The XBraid setup involves initializing the <kbd>app</kbd> struct and configuring XBraid for the desired number of timesteps, number of iterations, etc.
+The functions implemented for XBraid's use are declared in <kbd>BraidFuncs.hh</kbd> and defined in <kbd>BraidFuncs.cc</kbd>. The <kbd>HeatEquation</kbd> class and all <kbd>deal.ii</kbd> functionality is declared in <kbd>HeatEquation.hh</kbd> and defiend in <kbd>HeatEquationImplem.hh</kbd>. Since <kbd>HeatEquation</kbd> is a class template, its definition file <kbd>HeatEquationImplem.hh</kbd> is included at the bottom of <kbd>HeatEquation.hh</kbd>. Lastly various helper functions and variables such as the current processor id and the output stream are declared in <kbd>Utilities.hh</kbd> and defined in <kbd>Utilities.cc</kbd>.
+
+#### The test directory
+
+This directory contains tests to be built and run with CMake.
+These tests verify the correct implementation of the various functions.
+
+#### The doc directory
+
+This directory is for storing further documentation of the code.
+Not much is in this directory right now as most of the documentation is in this file (Readme.md) or in comments in the source code files.
+
+## Compiling
+
+To compile, you need deal.ii and XBraid to be installed with development headers somwehere on your system.
+Some implementation of MPI such as OpenMPI with development headers must also be installed.
+The source code for deal.ii is available at [https://dealii.org/](https://dealii.org/) and the source code for XBraid is available at [https://computation.llnl.gov/projects/parallel-time-integration-multigrid](https://computation.llnl.gov/projects/parallel-time-integration-multigrid).
+See the documentation of each package for compilation and installation instructions.
+
+Depending on where they are installed, parallel_in_time may need help finding these libraries.
+To find deal.ii, parallel_in_time first looks in typical deal.ii install directories followed by one directory up(<kbd>../</kbd>), two directories up (<kbd>../../</kbd>), and lastly in the environment variable <kbd>DEAL_II_DIR</kbd>.
+In contrast, XBraid currently does not have any default locations to look for and so the environment variable <kbd>BRAID_DIR</kbd> must be specified.
+For MPI, parallel_in_time looks in standard installation folders only, for that reason I recommend you install MPI with your package manager.
+
+A compile process of the parallel in time code may look like,
+
+ mkdir build
+ cd build
+ BRAID_DIR=/path/to/braid/ cmake ../
+ make
+
+There is currently no option to install parallel_in_time anywhere.
+The binaries are generated in the <kbd>bin</kbd> folder, and tests are placed into the <kbd>test</kbd> folder.
+Options that can be passed to CMake for parallel_in_time include:
+
+* <kbd>CMAKE_BUILD_TYPE=Debug/Release</kbd>
+* <kbd>DO_MFG=ON/OFF</kbd>
+* <kbd>USE_MPI=ON/OFF</kbd>
+
+The build type specifies whether to compile with debugging symbols, assertions, and optimizations or not.
+
+The option for manufactured solutions (<kbd>DO_MFG</kbd>) switches from solving the "standard" heat equation to solving a known solution heat equation so that the correctness of the code can be tested.
+
+Lastly the MPI option is only used to specify where to write output information when using the <kbd>pout()</kbd> function from the Utilities.hh file.
+If <kbd>USE_MPI</kbd> is set to <kbd>ON</kbd>, then every processor writes to its own file called <kbd>pout.<#></kbd> where <kbd><#></kbd> is the processor number.
+If <kbd>USE_MPI</kbd> is set to <kbd>OFF</kbd>, then every processor writes to stdout.
+
+## Running
+
+Once parallel_in_time has been compiled, the program can be run by calling the binary generated in <kbd>./build/bin/</kbd>.
+The test can be run by calling <kbd>ctest</kbd> from inside the <kbd>build/</kbd> directory.
+Unless the output path has been changed in the source code (currently hardcoded), then the output files will be placed into the folder the command was called from.
+
+There are no argument parameters for any of the binaries or tests.
+
+## Results
+
+To test the performance, a strong scaling study is run.
+The spatial grid is fixed at 3201 degrees of freedom, and the spatial grid consists of 25,000 uniform timesteps.
+No spatial parallelization is used and the grid is fixed for all timesteps.
+The parallel in time solution is solved using XBraid's multigrid reduction in time algorithm on 1, 2, 4, 16, 32, and 64 processors.
+The serial in time solution is run on a single processor using traditional sequential time stepping.
+The results are shown in the figure below.
+Running the multigrid algorithm on a single processor takes about an order of magnitude longer to run on a single processor than the serial algorihtm on a single processor.
+At 16 processors the multigrid algorithm the wall clock time is approximately the same for the serial algorithm as for the multigrid algorithm, and for 32 and 64 processors in time the wall clock is faster by about a factor of 2 for 64 processors.
+
+
+
+## Conclusions
+
+Using 64 times as many proceossors results in a speedup factor of approximately 2.
+This is a fairly heavy computational cost for a slight speedup.
+For comparison, in the reference paper Falgout et. al. they achieved a speedup of approximately 10 times when using 256 times as many processors as their serial case when solving the heat equation in two dimensions.
+A similar increase in processors may be needed to increase the speedup factor of this code from 2 to 10.
+The choice of whether to use serial in time or parallel in time largely rests in the size of problem being solved and the amount of computing resources available. Increasing the required number of timesteps will benefit the parallel in time algorithm provided enough extra resources are available.
+
+## Future Work
+
+There are many routes available for future work.
+First there are many optimizations that could be made to speed up the existing code base such as exploring the different multigrid cycles and finding the optimal load balancing.
+Ultimately those optimizations will probably only result in marginal gains per the Falgout paper.
+Allowing XBraid to prolong and restrict the spatial grid may be one of the more promising avenues of improvement.
+
+Future work that is of interest to the authors of XBraid is the development of adaptive mesh refinement (AMR) in a parallel in time algorithm.
+A particular challenges with parallel-in-time AMR is the time subcylcing that occurs in spatial subdomains with sequential time stepping.
+This code base does not use spatial subdomains with AMR and so could provide an easier understanding of the coarsening and refining of the time domain.
--- /dev/null
+Joshua Christopher <chrisjoshtopher@gmail.com>
--- /dev/null
+DEAL_II_WITH_MPI
--- /dev/null
+Parallel in Time Heat Equation
--- /dev/null
+Solves the heat equation time steps in parallel using the multigrid reduction in time technique from the XBraid library.
--- /dev/null
+/*-------- Project --------*/
+#include "BraidFuncs.hh"
+
+int my_Step(braid_App app,
+ braid_Vector ustop,
+ braid_Vector fstop,
+ braid_Vector u,
+ braid_StepStatus status)
+{
+ UNUSED(ustop);
+ UNUSED(fstop);
+ double tstart; /* current time */
+ double tstop; /* evolve to this time*/
+ int level;
+ double deltaT;
+
+ int index;
+ braid_StepStatusGetLevel(status, &level);
+ braid_StepStatusGetTstartTstop(status, &tstart, &tstop);
+ braid_StepStatusGetTIndex(status, &index);
+
+ deltaT = tstop - tstart;
+
+ dealii::Vector<double>& solution = u->data;
+
+ HeatEquation<2>& heateq = app->eq;
+
+ heateq.step(solution, deltaT, tstart, index);
+
+ return 0;
+}
+
+int
+my_Init(braid_App app,
+ double t,
+ braid_Vector *u_ptr)
+{
+ my_Vector *u = new(my_Vector);
+ int size = app->eq.size();
+ u->data.reinit(size);
+
+ app->eq.initialize(t, u->data);
+
+ *u_ptr = u;
+
+ return 0;
+}
+
+int
+my_Clone(braid_App app,
+ braid_Vector u,
+ braid_Vector *v_ptr)
+{
+ UNUSED(app);
+ my_Vector *v = new(my_Vector);
+ int size = u->data.size();
+ v->data.reinit(size);
+ for(size_t i=0, end=v->data.size(); i != end; ++i)
+ {
+ v->data[i] = u->data[i];
+ }
+ *v_ptr = v;
+
+ return 0;
+}
+
+int
+my_Free(braid_App app,
+ braid_Vector u)
+{
+ UNUSED(app);
+ delete u;
+
+ return 0;
+}
+
+int my_Sum(braid_App app,
+ double alpha,
+ braid_Vector x,
+ double beta,
+ braid_Vector y)
+{
+ UNUSED(app);
+ Vector<double>& vec = y->data;
+ vec.sadd(beta, alpha, x->data);
+
+ return 0;
+}
+
+int
+my_SpatialNorm(braid_App app,
+ braid_Vector u,
+ double *norm_ptr)
+{
+ UNUSED(app);
+ double dot = 0.0;
+ dot = u->data.l2_norm();
+ *norm_ptr = dot;
+
+ return 0;
+}
+
+int
+my_Access(braid_App app,
+ braid_Vector u,
+ braid_AccessStatus astatus)
+{
+ double t;
+ int index;
+
+ braid_AccessStatusGetT(astatus, &t);
+ braid_AccessStatusGetTIndex(astatus, &index);
+
+ app->eq.output_results(index, t, u->data);
+
+#if DO_MFG
+ if(index == app->final_step)
+ {
+ pout() << "Doing error calc of step: " << index << std::endl;
+ app->eq.process_solution(t, index, u->data);
+ }
+#endif
+
+ return 0;
+}
+
+int
+my_BufSize(braid_App app,
+ int *size_ptr,
+ braid_BufferStatus bstatus)
+{
+ UNUSED(bstatus);
+ int size = app->eq.size();
+ *size_ptr = (size+1)*sizeof(double);
+
+ return 0;
+}
+
+int
+my_BufPack(braid_App app,
+ braid_Vector u,
+ void *buffer,
+ braid_BufferStatus bstatus)
+{
+
+ UNUSED(app);
+ double *dbuffer = (double*)buffer;
+ int size = u->data.size();
+ dbuffer[0] = size;
+ for(int i=0; i != size; ++i)
+ {
+ dbuffer[i+1] = (u->data)[i];
+ }
+ braid_BufferStatusSetSize(bstatus, (size+1)*sizeof(double));
+
+ return 0;
+}
+
+int
+my_BufUnpack(braid_App app,
+ void *buffer,
+ braid_Vector *u_ptr,
+ braid_BufferStatus bstatus)
+{
+ UNUSED(app);
+ UNUSED(bstatus);
+
+ my_Vector *u = NULL;
+ double *dbuffer = (double*)buffer;
+ int size = dbuffer[0];
+ u = new(my_Vector);
+ u->data.reinit(size);
+
+ for(int i = 0; i != size; ++i)
+ {
+ (u->data)[i] = dbuffer[i+1];
+ }
+ *u_ptr = u;
+
+ return 0;
+}
--- /dev/null
+#ifndef _BRAIDFUNCS_H_
+#define _BRAIDFUNCS_H_
+
+/**
+ * \file BraidFuncs.cc
+ * \brief Contains the implementation of the mandatory X-Braid functions
+ *
+ * X-Braid mandates several functions in order to drive the solution.
+ * This file contains the implementation of said mandatory functions.
+ * See the X-Braid documentation for more information.
+ * There are several functions that are optional in X-Braid that may
+ * or may not be implemented in here.
+ *
+ */
+
+/*-------- System --------*/
+
+/*-------- Third Party --------*/
+#include <deal.II/numerics/vector_tools.h>
+
+#include "braid.h"
+#include "braid_test.h"
+
+/*-------- Project --------*/
+#include "HeatEquation.hh"
+
+/**
+ * \brief Struct that contains the data
+ */
+typedef struct _braid_Vector_struct
+{
+ dealii::Vector<double> data;
+} my_Vector;
+
+// Wrap the Heat Equation in a struct
+typedef struct _braid_App_struct
+{
+ HeatEquation<2> eq;
+ int final_step;
+} my_App;
+
+
+/**
+ * @brief my_Step - Takes a step in time, advancing the u vector
+ *
+ * @param app - The braid app struct
+ * @param ustop - The solution data at the end of this time step
+ * @param fstop - RHS data (such as forcing function?)
+ * @param u - The solution data at the beginning of this time step
+ * @param status - Status structure that contains various info of this time
+ *
+ * @return Success (0) or failure (1)
+ **/
+int my_Step(braid_App app,
+ braid_Vector ustop,
+ braid_Vector fstop,
+ braid_Vector u,
+ braid_StepStatus status);
+
+
+/**
+ * @brief my_Init - Initializes a solution data at the given time
+ * For now, initializes the solution to zero no matter what time we are at
+ *
+ * @param app - The braid app struct containing user data
+ * @param t - Time at which the solution is initialized
+ * @param u_ptr - The solution data that needs to be filled
+ *
+ * @return Success (0) or failure (1)
+ **/
+int
+my_Init(braid_App app,
+ double t,
+ braid_Vector *u_ptr);
+
+
+/**
+ * @brief my_Clone - Clones a vector into a new vector
+ *
+ * @param app - The braid app struct containing user data
+ * @param u - The existing vector containing data
+ * @param v_ptr - The empty vector that needs to be filled
+ *
+ * @return Success (0) or failure (1)
+ **/
+int
+my_Clone(braid_App app,
+ braid_Vector u,
+ braid_Vector *v_ptr);
+
+
+/**
+ * @brief my_Free - Deletes a vector
+ *
+ * @param app - The braid app struct containing user data
+ * @param u - The vector that needs to be deleted
+ *
+ * @return Success (0) or failure (1)
+ **/
+int
+my_Free(braid_App app,
+ braid_Vector u);
+
+
+/**
+ * @brief my_Sum - Sums two vectors in an AXPY operation
+ * The operation is y = alpha*x + beta*y
+ *
+ * @param app - The braid app struct containing user data
+ * @param alpha - The coefficient in front of x
+ * @param x - A vector that is multiplied by alpha then added to y
+ * @param beta - The coefficient of y
+ * @param y - A vector that is multiplied by beta then summed with x
+ *
+ * @return Success (0) or failure (1)
+ **/
+int
+my_Sum(braid_App app,
+ double alpha,
+ braid_Vector x,
+ double beta,
+ braid_Vector y);
+
+/**
+ * \brief Returns the spatial norm of the provided vector
+ *
+ * Calculates and returns the spatial norm of the provided vector.
+ * Interestingly enough, X-Braid does not specify a particular norm.
+ * to keep things simple, we implement the Euclidean norm.
+ *
+ * \param app - The braid app struct containing user data
+ * \param u - The vector we need to take the norm of
+ * \param norm_ptr - Pointer to the norm that was calculated, need to modify this
+ * \return Success (0) or failure (1)
+ */
+int
+my_SpatialNorm(braid_App app,
+ braid_Vector u,
+ double *norm_ptr);
+
+/**
+ * \brief Allows the user to output details
+ *
+ * The Access function is called at various points to allow the user to output
+ * information to the screen or to files.
+ * The astatus parameter provides various information about the simulation,
+ * see the XBraid documentation for details on what information you can get.
+ * Example information is what the current timestep number and current time is.
+ * If the access level (in parallel_in_time.cc) is set to 0, this function is
+ * never called.
+ * If the access level is set to 1, the function is called after the last
+ * XBraid cycle.
+ * If the access level is set to 2, it is called every XBraid cycle.
+ *
+ * \param app - The braid app struct containing user data
+ * \param u - The vector containing the data at the status provided
+ * \param astatus - The Braid status structure
+ * \return Success (0) or failure (1)
+ */
+int
+my_Access(braid_App app,
+ braid_Vector u,
+ braid_AccessStatus astatus);
+
+/**
+ * \brief Calculates the size of a buffer for MPI data transfer
+ *
+ * Calculates the size of the buffer that is needed to transfer
+ * a solution vector to another processor.
+ * The bstatus parameter provides various information on the
+ * simulation, see the XBraid documentation for all possible
+ * fields.
+ *
+ * \param app - The braid app struct containing user data
+ * \param size_ptr A pointer to the calculated size
+ * \param bstatus The XBraid status structure
+ * \return Success (0) or failure (1)
+ */
+int
+my_BufSize(braid_App app,
+ int *size_ptr,
+ braid_BufferStatus bstatus);
+
+/**
+ * \brief Linearizes a vector to be sent to another processor
+ *
+ * Linearizes (packs) a data buffer with the contents of
+ * some solution state u.
+ *
+ * \param app - The braid app struct containing user data
+ * \param u The vector that must be packed into buffer
+ * \param buffer The buffer that must be filled with u
+ * \param bstatus The XBraid status structure
+ * \return Success (0) or failure (1)
+ */
+int
+my_BufPack(braid_App app,
+ braid_Vector u,
+ void *buffer,
+ braid_BufferStatus bstatus);
+
+/**
+ * \brief Unpacks a vector that was sent from another processor
+ *
+ * Unpacks a linear data buffer into the vector pointed to by
+ * u_ptr.
+ *
+ * \param app - The braid app struct containing user data
+ * \param buffer The buffer that must be unpacked
+ * \param u_ptr The pointer to the vector that is filled
+ * \param bstatus The XBraid status structure
+ * \return Success (0) or failure (1)
+ */
+int
+my_BufUnpack(braid_App app,
+ void *buffer,
+ braid_Vector *u_ptr,
+ braid_BufferStatus bstatus);
+
+#endif // _BRAIDFUNCS_H_
--- /dev/null
+#ifndef _HEATEQUATION_H_
+#define _HEATEQUATION_H_
+
+#include <fstream>
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/base/convergence_table.h>
+
+using namespace dealii;
+
+template <int dim>
+class HeatEquation
+{
+public:
+ HeatEquation();
+ void define();
+ void step(Vector<double>& braid_data,
+ double deltaT,
+ double a_time,
+ int a_time_idx);
+
+ int size() const; /// Returns the size of the solution vector
+
+ void dump_vec(std::ofstream& file,
+ const Vector<double>& vector) const;
+
+ void output_results(int a_time_idx,
+ double a_time,
+ Vector<double>& a_solution) const;
+
+ void initialize(double a_time,
+ Vector<double>& a_vector) const;
+
+ void process_solution(double a_time,
+ int a_index,
+ const Vector<double>& a_vector);
+
+private:
+ void setup_system();
+ void solve_time_step(Vector<double>& a_solution);
+
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> mass_matrix;
+ SparseMatrix<double> laplace_matrix;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> system_rhs;
+
+ std::ofstream myfile;
+
+ const double theta;
+
+ // These were originally in the run() function but because
+ // I am splitting the run() function up into define and step
+ // they need to become member data
+ Vector<double> tmp;
+ Vector<double> forcing_terms;
+
+ ConvergenceTable convergence_table;
+};
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ RightHandSide ()
+ :
+ Function<dim>(),
+ period (0.2)
+ {}
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+private:
+ const double period;
+};
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+template <int dim>
+class RightHandSideMFG : public Function<dim>
+{
+public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+template <int dim>
+class InitialValuesMFG : public Function<dim>
+{
+public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+/**
+ * Provides the exact value for the manufactured solution.
+ * User must set the time using ExactValuesMFG::set_time(const Number new_time);
+ * Calculates \f$u_h=e^(-4\pi\pi t)*\cos(2\pi x)\cos(2\pi y)\f$ for all points given to it.
+ * In the context of manufactured solutions, this is used on the boundary.
+ */
+template <int dim>
+class ExactValuesMFG : public Function<dim>
+{
+public:
+ /**
+ * \brief Computes the value at the given point and member data time
+ *
+ * Computes the exact value of the manufactured solution at point p and
+ * the member data time. See the class documentation and the design doc
+ * for details on what the exact solution is.
+ *
+ * \param p The point that the exact solution is computed at
+ * \param component The component of the exact solution (always 0 for now)
+ * \return double The exact value that was computed
+ */
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ /**
+ * \brief Computes the gradient of the exact solution at the given point
+ *
+ * Computes the gradient of the exact/manufactured solution value at
+ * point p and member data time. See the design doc for details on
+ * what the gradient of the exact solution is
+ *
+ * \param p The point that the gradient is calculated at
+ * \param component The component of the system of equations this gradient is for
+ * \return Tensor<1,dim> A rank 1 tensor that contains the gradient
+ * in each spatial dimension
+ */
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+#include "HeatEquationImplem.hh"
+
+#endif // _HEATEQUATION_H_
--- /dev/null
+#include <iomanip>
+#include <math.h>
+#include "Utilities.hh"
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ (void) component;
+ Assert (component == 0, ExcIndexRange(component, 0, 1));
+ Assert (dim == 2, ExcNotImplemented());
+
+ double time = this->get_time();
+
+ if ((p[0] > 0.5) && (p[1] > -0.5))
+ {
+ return std::exp(-0.5*(time-0.125)*(time-0.125)/(0.005));
+ }
+ else if ((p[0] > -0.5) && (p[1] > 0.5))
+ {
+ return std::exp(-0.5*(time-0.375)*(time-0.375)/(0.005));
+ }
+ else
+ {
+ return 0;
+ }
+
+ return 0; // No forcing function
+}
+
+template <int dim>
+double RightHandSideMFG<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ (void) component;
+ Assert (component == 0, ExcIndexRange(component, 0, 1));
+ Assert (dim == 2, ExcNotImplemented());
+
+ double time = this->get_time();
+
+ // return the manufactured solution of the right hand side
+ double pi = numbers::PI;
+ return 4*pi*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
+}
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ UNUSED(p);
+ (void) component;
+ Assert (component == 0, ExcIndexRange(component, 0, 1));
+ return 0;
+}
+
+template <int dim>
+double ExactValuesMFG<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ (void) component;
+ Assert (component == 0, ExcIndexRange(component, 0, 1));
+
+ double time = this->get_time();
+ const double pi = numbers::PI;
+ // Return our manufactured solution boundary value
+ return std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
+}
+
+
+// TODO: Find bug in here
+template <int dim>
+Tensor<1,dim> ExactValuesMFG<dim>::gradient (const Point<dim> &p,
+ const unsigned int) const
+{
+ Tensor<1,dim> return_value;
+ const double pi = numbers::PI;
+ double time = this->get_time();
+ return_value[0] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[1])*std::sin(2*pi*p[0]);
+ return_value[1] = -2*pi*std::exp(-4*pi*pi*time)*std::cos(2*pi*p[0])*std::sin(2*pi*p[1]);
+ return return_value;
+}
+
+template <int dim>
+double InitialValuesMFG<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ (void) component;
+ Assert (component == 0, ExcIndexRange(component, 0, 1));
+ const double pi = numbers::PI;
+ // Return our manufactured solution initial value
+ return std::cos(2*pi*p[0])*std::cos(2*pi*p[1]);
+}
+
+template <int dim>
+HeatEquation<dim>::HeatEquation ()
+ :
+ fe(1),
+ dof_handler(triangulation),
+ theta(0.5)
+{
+}
+
+template <int dim>
+void HeatEquation<dim>::initialize(double a_time,
+ Vector<double>& a_vector) const
+{
+#if DO_MFG
+ // We only initialize values in the manufactured solution case
+ InitialValuesMFG<dim> iv_function;
+ iv_function.set_time(a_time);
+ VectorTools::project (dof_handler, constraints,
+ QGauss<dim>(fe.degree+1), iv_function,
+ a_vector);
+#else
+ UNUSED(a_time);
+ UNUSED(a_vector);
+#endif // DO_MFG
+ // If not the MFG solution case, a_vector is already zero'd so do nothing
+}
+
+template <int dim>
+void HeatEquation<dim>::setup_system()
+{
+ dof_handler.distribute_dofs(fe);
+
+ // pout() << std::endl
+ // << "==========================================="
+ // << std::endl
+ // << "Number of active cells: " << triangulation.n_active_cells()
+ // << std::endl
+ // << "Number of degrees of freedom: " << dof_handler.n_dofs()
+ // << std::endl
+ // << std::endl;
+
+ constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+ constraints.close();
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler,
+ dsp,
+ constraints,
+ /*keep_constrained_dofs = */ true);
+ sparsity_pattern.copy_from(dsp);
+
+ mass_matrix.reinit(sparsity_pattern);
+ laplace_matrix.reinit(sparsity_pattern);
+ system_matrix.reinit(sparsity_pattern);
+
+ MatrixCreator::create_mass_matrix(dof_handler,
+ QGauss<dim>(fe.degree+1),
+ mass_matrix);
+ MatrixCreator::create_laplace_matrix(dof_handler,
+ QGauss<dim>(fe.degree+1),
+ laplace_matrix);
+
+ system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void HeatEquation<dim>::solve_time_step(Vector<double>& a_solution)
+{
+ SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
+ SolverCG<> cg(solver_control);
+
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.0);
+
+ cg.solve(system_matrix, a_solution, system_rhs,
+ preconditioner);
+
+ constraints.distribute(a_solution);
+}
+
+
+
+template <int dim>
+void HeatEquation<dim>::output_results(int a_time_idx,
+ double a_time,
+ Vector<double>& a_solution) const
+{
+
+ DataOutBase::VtkFlags vtk_flags;
+ vtk_flags.time = a_time;
+ vtk_flags.cycle = a_time_idx;
+
+ DataOut<dim> data_out;
+ data_out.set_flags(vtk_flags);
+
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(a_solution, "U");
+
+ data_out.build_patches();
+
+ const std::string filename = "solution-"
+ + Utilities::int_to_string(a_time_idx, 3) +
+ ".vtk";
+ std::ofstream output(filename.c_str());
+ data_out.write_vtk(output);
+}
+
+// This function won't make much sense in real parallel in time...
+template <int dim>
+void HeatEquation<dim>::define()
+{
+ const unsigned int initial_global_refinement = 6;
+
+ GridGenerator::hyper_L (triangulation);
+ triangulation.refine_global (initial_global_refinement);
+
+ setup_system();
+
+ tmp.reinit (dof_handler.n_dofs());
+ forcing_terms.reinit (dof_handler.n_dofs());
+}
+
+template<int dim>
+void HeatEquation<dim>::step(Vector<double>& braid_data,
+ double deltaT,
+ double a_time,
+ int a_time_idx)
+{
+ // Set old solution to the braid data
+ // old_solution = braid_data;
+ a_time += deltaT;
+ ++a_time_idx;
+
+ mass_matrix.vmult(system_rhs, braid_data);
+
+ laplace_matrix.vmult(tmp, braid_data);
+
+ system_rhs.add(-(1 - theta) * deltaT, tmp);
+
+#if DO_MFG
+ RightHandSideMFG<dim> rhs_function;
+#else
+ RightHandSide<dim> rhs_function;
+#endif
+ rhs_function.set_time(a_time);
+ VectorTools::create_right_hand_side(dof_handler,
+ QGauss<dim>(fe.degree+1),
+ rhs_function,
+ tmp);
+
+ forcing_terms = tmp;
+ forcing_terms *= deltaT * theta;
+
+ rhs_function.set_time(a_time - deltaT);
+ VectorTools::create_right_hand_side(dof_handler,
+ QGauss<dim>(fe.degree+1),
+ rhs_function,
+ tmp);
+
+ forcing_terms.add(deltaT * (1 - theta), tmp);
+ system_rhs += forcing_terms;
+
+ system_matrix.copy_from(mass_matrix);
+ system_matrix.add(theta * deltaT, laplace_matrix);
+
+ constraints.condense (system_matrix, system_rhs);
+
+ {
+#if DO_MFG
+ // Set boundary to exact value in MFG solution
+ ExactValuesMFG<dim> boundary_values_function;
+#else
+ BoundaryValues<dim> boundary_values_function;
+#endif
+ boundary_values_function.set_time(a_time);
+
+ std::map<types::global_dof_index, double> boundary_values;
+ VectorTools::interpolate_boundary_values(dof_handler,
+ 0,
+ boundary_values_function,
+ boundary_values);
+
+ MatrixTools::apply_boundary_values(boundary_values,
+ system_matrix,
+ braid_data,
+ system_rhs);
+ }
+
+ solve_time_step(braid_data);
+}
+
+template<int dim>
+int HeatEquation<dim>::size() const
+{
+ return dof_handler.n_dofs();
+}
+
+template<int dim> void
+HeatEquation<dim>::dump_vec(std::ofstream& file,
+ const Vector<double>& vector) const
+{
+ file << "Dumping vec:";
+ for(int i=0; i != vector.size(); ++i)
+ {
+ file << "\n" << vector[i];
+ }
+ file << std::endl;
+}
+
+template<int dim> void
+HeatEquation<dim>::process_solution(double a_time,
+ int a_index,
+ const Vector<double>& a_vector)
+{
+ // Compute the exact value for the manufactured solution case
+ ExactValuesMFG<dim> exact_function;
+ exact_function.set_time(a_time);
+
+ Vector<double> difference_per_cell (triangulation.n_active_cells());
+ VectorTools::integrate_difference(dof_handler,
+ a_vector,
+ exact_function,
+ difference_per_cell,
+ QGauss<dim>(fe.degree+1),
+ VectorTools::L2_norm);
+
+ const double L2_error = VectorTools::compute_global_error(triangulation,
+ difference_per_cell,
+ VectorTools::L2_norm);
+
+ VectorTools::integrate_difference(dof_handler,
+ a_vector,
+ exact_function,
+ difference_per_cell,
+ QGauss<dim>(fe.degree+1),
+ VectorTools::H1_seminorm);
+
+ const double H1_error = VectorTools::compute_global_error(triangulation,
+ difference_per_cell,
+ VectorTools::H1_seminorm);
+
+ const QTrapez<1> q_trapez;
+ const QIterated<dim> q_iterated (q_trapez, 5);
+ VectorTools::integrate_difference (dof_handler,
+ a_vector,
+ exact_function,
+ difference_per_cell,
+ q_iterated,
+ VectorTools::Linfty_norm);
+ const double Linfty_error = VectorTools::compute_global_error(triangulation,
+ difference_per_cell,
+ VectorTools::Linfty_norm);
+
+ const unsigned int n_active_cells = triangulation.n_active_cells();
+ const unsigned int n_dofs = dof_handler.n_dofs();
+
+ std::cout << "Cycle " << a_index << ':'
+ << std::endl
+ << " Number of active cells: "
+ << n_active_cells
+ << std::endl
+ << " Number of degrees of freedom: "
+ << n_dofs
+ << std::endl;
+
+ convergence_table.add_value("cycle", a_index);
+ convergence_table.add_value("cells", n_active_cells);
+ convergence_table.add_value("dofs", n_dofs);
+ convergence_table.add_value("L2", L2_error);
+ convergence_table.add_value("H1", H1_error);
+ convergence_table.add_value("Linfty", Linfty_error);
+
+ convergence_table.set_precision("L2", 3);
+ convergence_table.set_precision("H1", 3);
+ convergence_table.set_precision("Linfty", 3);
+
+ convergence_table.set_scientific("L2", true);
+ convergence_table.set_scientific("H1", true);
+ convergence_table.set_scientific("Linfty", true);
+
+ convergence_table.set_tex_caption("cells", "\\# cells");
+ convergence_table.set_tex_caption("dofs", "\\# dofs");
+ convergence_table.set_tex_caption("L2", "@f$L^2@f$-error");
+ convergence_table.set_tex_caption("H1", "@f$H^1@f$-error");
+ convergence_table.set_tex_caption("Linfty", "@f$L^\\infty@f$-error");
+
+ convergence_table.set_tex_format("cells", "r");
+ convergence_table.set_tex_format("dofs", "r");
+
+ std::cout << std::endl;
+ convergence_table.write_text(std::cout);
+
+ std::ofstream error_table_file("tex-conv-table.tex");
+ convergence_table.write_tex(error_table_file);
+}
--- /dev/null
+#include "Utilities.hh"
+
+#include <string>
+#include <fstream>
+
+#include "mpi.h"
+
+int procID = 0;
+
+// the shared variables
+
+static std::string s_pout_filename ;
+static std::string s_pout_basename ;
+static std::ofstream s_pout ;
+
+static bool s_pout_init = false ;
+static bool s_pout_open = false ;
+
+#ifdef USE_MPI
+// in parallel, compute the filename give the basename
+//[NOTE: dont call this before MPI is initialized.]
+static void setFileName()
+{
+ static const size_t ProcnumSize = 1 + 10 + 1 ; //'.' + 10digits + '\0'
+ char procnum[ProcnumSize] ;
+ snprintf( procnum ,ProcnumSize ,".%d" ,procID);
+ s_pout_filename = s_pout_basename + procnum ;
+}
+
+// in parallel, close the file if nec., open it and check for success
+static void openFile()
+{
+ if ( s_pout_open )
+ {
+ s_pout.close();
+ }
+ s_pout.open( s_pout_filename.c_str() );
+ // if open() fails, we have problems, but it's better
+ // to try again later than to make believe it succeeded
+ s_pout_open = (bool)s_pout ;
+}
+
+#else
+// in serial, filename is always cout
+static void setFileName()
+{
+ s_pout_filename = "cout" ;
+}
+
+// in serial, this does absolutely nothing
+static void openFile()
+{
+}
+#endif
+
+std::ostream& pout()
+{
+#ifdef USE_MPI
+ // the common case is _open == true, which just returns s_pout
+ if ( ! s_pout_open )
+ {
+ // the uncommon cae: the file isn't opened, MPI may not be
+ // initialized, and the basename may not have been set
+ int flag_i, flag_f;
+ MPI_Initialized(&flag_i);
+ MPI_Finalized(&flag_f);
+ // app hasn't set a basename yet, so set the default
+ if ( ! s_pout_init )
+ {
+ s_pout_basename = "pout" ;
+ s_pout_init = true ;
+ }
+ // if MPI not initialized, we cant open the file so return cout
+ if ( ! flag_i || flag_f)
+ {
+ return std::cout; // MPI hasn't been started yet, or has ended....
+ }
+ // MPI is initialized, so file must not be, so open it
+ setFileName() ;
+ openFile() ;
+ // finally, in case the open failed, return cout
+ if ( ! s_pout_open )
+ {
+ return std::cout ;
+ }
+ }
+ return s_pout ;
+#else
+ return std::cout;
+#endif
+}
--- /dev/null
+#ifndef _UTILITIES_H_
+#define _UTILITIES_H_
+#include <iostream>
+
+#define UNUSED(x) (void)(x)
+
+extern int procID;
+
+std::ostream& pout();
+#endif
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2013 - 2018 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+
+ *
+ * Author: Joshua Christopher, Colorado State University, 2018
+ */
+
+#include <fstream>
+#include <iostream>
+
+#include "BraidFuncs.hh"
+#include "HeatEquation.hh"
+#include "Utilities.hh"
+
+int main(int argc, char *argv[])
+{
+ try
+ {
+ using namespace dealii;
+
+ /* Initialize MPI */
+ MPI_Comm comm; //, comm_x, comm_t;
+ int rank;
+ MPI_Init(&argc, &argv);
+ comm = MPI_COMM_WORLD;
+ MPI_Comm_rank(comm, &rank);
+ procID = rank;
+
+ // Set up X-Braid
+ /* Initialize Braid */
+ braid_Core core;
+ double tstart = 0.0;
+ double tstop = 0.002;
+ double ntime = 10;
+ my_App *app = new(my_App);
+
+ braid_Init(MPI_COMM_WORLD, comm, tstart, tstop, ntime, app,
+ my_Step, my_Init, my_Clone, my_Free, my_Sum, my_SpatialNorm,
+ my_Access, my_BufSize, my_BufPack, my_BufUnpack, &core);
+
+ /* Define XBraid parameters
+ * See -help message forf descriptions */
+ int max_levels = 3;
+ // int nrelax = 1;
+ // int skip = 0;
+ double tol = 1.e-7;
+ // int cfactor = 2;
+ int max_iter = 5;
+ // int min_coarse = 10;
+ // int fmg = 0;
+ // int scoarsen = 0;
+ // int res = 0;
+ // int wrapper_tests = 0;
+ int print_level = 1;
+ int access_level = 1;
+ int use_sequential= 0;
+
+ braid_SetPrintLevel( core, print_level);
+ braid_SetAccessLevel( core, access_level);
+ braid_SetMaxLevels(core, max_levels);
+ // braid_SetMinCoarse( core, min_coarse );
+ // braid_SetSkip(core, skip);
+ // braid_SetNRelax(core, -1, nrelax);
+ braid_SetAbsTol(core, tol);
+ // braid_SetCFactor(core, -1, cfactor);
+ braid_SetMaxIter(core, max_iter);
+ braid_SetSeqSoln(core, use_sequential);
+
+ app->eq.define();
+ app->final_step = ntime;
+
+ braid_Drive(core);
+
+ // Free the memory now that we are done
+ braid_Destroy(core);
+
+ delete app;
+
+ // Clean up MPI
+ // MPI_Comm_free(&comm);
+ MPI_Finalize();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl << exc.what()
+ << std::endl << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl << "Aborting!"
+ << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
+
--- /dev/null
+set(TEST_LIBRARIES ${LIB_NAME} ${BRAID_DIR}/libbraid.a)
+DEAL_II_PICKUP_TESTS()
--- /dev/null
+#include <iostream>
+
+#include "braid.h"
+#include "braid_test.h"
+
+#include "BraidFuncs.hh"
+
+int main(int argc, char** argv)
+{
+ MPI_Comm comm;
+ int rank;
+ MPI_Init(&argc, &argv);
+ comm = MPI_COMM_WORLD;
+ MPI_Comm_rank(comm, &rank);
+
+ my_App *app = new(my_App);
+ app->eq.define();
+
+ double time = 0.2;
+
+ braid_Int init_access_result = braid_TestInitAccess(app,
+ comm,
+ stdout,
+ time,
+ my_Init,
+ my_Access,
+ my_Free);
+
+ braid_Int clone_result = braid_TestClone(app,
+ comm,
+ stdout,
+ time,
+ my_Init,
+ my_Access,
+ my_Free,
+ my_Clone);
+
+ braid_Int sum_result = braid_TestSum(app,
+ comm,
+ stdout,
+ time,
+ my_Init,
+ my_Access,
+ my_Free,
+ my_Clone,
+ my_Sum);
+
+ braid_Int norm_result = braid_TestSpatialNorm(app,
+ comm,
+ stdout,
+ time,
+ my_Init,
+ my_Free,
+ my_Clone,
+ my_Sum,
+ my_SpatialNorm);
+
+ braid_Int buf_result = braid_TestBuf(app,
+ comm,
+ stdout,
+ time,
+ my_Init,
+ my_Free,
+ my_Sum,
+ my_SpatialNorm,
+ my_BufSize,
+ my_BufPack,
+ my_BufUnpack);
+ // /* Create spatial communicator for wrapper-tests */
+ // braid_SplitCommworld(&comm, 1, &comm_x, &comm_t);
+//
+// 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);
+
+ /* Finalize MPI */
+ MPI_Finalize();
+}
--- /dev/null
+
+Starting braid_TestInitAccess
+
+ braid_TestInitAccess: Starting Test 1
+ braid_TestInitAccess: u = init(t=2.00e-01)
+ braid_TestInitAccess: access(u)
+ braid_TestInitAccess: check output: wrote u for initial condition at t=2.00e-01.
+
+ braid_TestInitAccess: free(u)
+Finished braid_TestInitAccess
+
+Starting braid_TestClone
+
+ braid_TestClone: Starting Test 1
+ braid_TestClone: u = init(t=2.00e-01)
+ braid_TestClone: v = clone(u)
+ braid_TestClone: access(u)
+ braid_TestClone: access(v)
+ braid_TestClone: check output: wrote u and v for initial condition at t=2.00e-01.
+
+ braid_TestClone: free(u)
+ braid_TestClone: free(v)
+Finished braid_TestClone
+
+Starting braid_TestSum
+
+ braid_TestSum: Starting Test 1
+ braid_TestSum: u = init(t=2.00e-01)
+ braid_TestSum: v = clone(u)
+ braid_TestSum: v = u - v
+ braid_TestSum: access(v)
+ braid_TestSum: check output: v should equal the zero vector
+
+ braid_TestSum: Starting Test 2
+ braid_TestSum: v = 2*u + v
+ braid_TestSum: access(v)
+ braid_TestSum: access(u)
+ braid_TestSum: check output: v should equal 2*u
+
+ braid_TestSum: free(u)
+ braid_TestSum: free(v)
+Finished braid_TestSum
+
+Starting braid_TestSpatialNorm
+
+ braid_TestSpatialNorm: Starting Test 1
+ braid_TestSpatialNorm: u = init(t=2.00e-01)
+ braid_TestSpatialNorm: spatialnorm(u)
+ braid_TestSpatialNorm: Warning: spatialnorm(u) = 0.0
+ braid_TestSpatialNorm: v = clone(u)
+ braid_TestSpatialNorm: v = u - v
+ braid_TestSpatialNorm: spatialnorm(v)
+ braid_TestSpatialNorm: Test 1 Passed
+ braid_TestSpatialNorm: actual output: spatialnorm(v) = 0.00e+00
+ braid_TestSpatialNorm: expected output: spatialnorm(v) = 0.0
+
+ braid_TestSpatialNorm: Starting Test 2
+ braid_TestSpatialNorm: w = clone(u)
+ braid_TestSpatialNorm: w = u + w
+ braid_TestSpatialNorm: spatialnorm(u)
+ braid_TestSpatialNorm: spatialnorm(w)
+ braid_TestSpatialNorm: Test 2 Failed, Likely due to u = 0
+ braid_TestSpatialNorm: actual output: spatialnorm(w) / spatialnorm(u) = 0.00e+00 / 0.00e+00 = -nan
+ braid_TestSpatialNorm: expected output: spatialnorm(w) / spatialnorm(u) = 2.0
+
+ braid_TestSpatialNorm: Starting Test 3
+ braid_TestSpatialNorm: free(w)
+ braid_TestSpatialNorm: w = clone(u)
+ braid_TestSpatialNorm: w = 0.0*u + 0.5*w
+ braid_TestSpatialNorm: spatialnorm(u)
+ braid_TestSpatialNorm: spatialnorm(w)
+ braid_TestSpatialNorm: Test 3 Failed, Likely due to u = 0
+ braid_TestSpatialNorm: actual output: spatialnorm(w) / spatialnorm(u) = 0.00e+00 / 0.00e+00 = -nan
+ braid_TestSpatialNorm: expected output: spatialnorm(w) / spatialnorm(u) = 0.5
+
+ braid_TestSpatialNorm: free(u)
+ braid_TestSpatialNorm: free(v)
+ braid_TestSpatialNorm: free(w)
+Finished braid_TestSpatialNorm: some tests failed, possibly due to u = 0
+
+Starting braid_TestBuf
+
+ braid_TestBuf: Starting Test 1
+ braid_TestBuf: u = init(t=2.00e-01)
+ braid_TestBuf: spatialnorm(u)
+ braid_TestBuf: Warning: spatialnorm(u) = 0.0
+ braid_TestBuf: size = bufsize()
+ braid_TestBuf: buffer = malloc(size)
+ braid_TestBuf: buffer = bufpack(u, buffer))
+ braid_TestBuf: v = bufunpack(buffer)
+ braid_TestBuf: v = u - v
+ braid_TestBuf: spatialnorm(v)
+ braid_TestBuf: Test 1 Passed
+ braid_TestBuf: actual output: spatialnorm(v) = 0.00e+00
+ braid_TestBuf: expected output: spatialnorm(v) = 0.0
+
+ braid_TestBuf: free(u)
+ braid_TestBuf: free(v)
+Finished braid_TestBuf: all tests passed successfully
projects / code-gallery.git / commitdiff
