diff --git a/.gitignore b/.gitignore
index 5dde0ae23..4ca2bb160 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,10 +2,8 @@ CMakeCache.txt
 CMakeFiles
 CMakeScripts
 *.cmake
-
-build*
+build/*
 TAGS
-
 Testing
 Makefile
 cmake_install.cmake
diff --git a/src/dg/dg.cpp b/src/dg/dg.cpp
index 7d71c7029..62fbbbd0c 100644
--- a/src/dg/dg.cpp
+++ b/src/dg/dg.cpp
@@ -55,7 +55,6 @@ ::create_discontinuous_galerkin(
     //} else if (pde_type == PDE_enum::convection_diffusion) {
     //    return new DG<dim,1,real>(parameters_input, degree);
     //}
-
     if (parameters_input->use_weak_form) {
         if (pde_type == PDE_enum::advection) {
             return std::make_shared< DGWeak<dim,1,real> >(parameters_input, degree);
@@ -91,15 +90,15 @@ ::create_discontinuous_galerkin(
 
 // DGBase ***************************************************************************
 template <int dim, typename real>
-DGBase<dim,real>::DGBase(
+DGBase<dim,real>::DGBase( // @suppress("Class members should be properly initialized")
     const int nstate_input,
     const Parameters::AllParameters *const parameters_input,
     const unsigned int max_degree_input)
-    : DGBase<dim,real>(nstate_input, parameters_input, max_degree_input, this->create_collection_tuple(max_degree_input, nstate_input))
+    : DGBase<dim,real>(nstate_input, parameters_input, max_degree_input, this->create_collection_tuple(max_degree_input, nstate_input, parameters_input))
 { }
 
 template <int dim, typename real>
-DGBase<dim,real>::DGBase(
+DGBase<dim,real>::DGBase( // @suppress("Class members should be properly initialized")
     const int nstate_input,
     const Parameters::AllParameters *const parameters_input,
     const unsigned int max_degree_input,
@@ -128,7 +127,7 @@ template <int dim, typename real>
 std::tuple< dealii::hp::MappingCollection<dim>, dealii::hp::FECollection<dim>,
             dealii::hp::QCollection<dim>, dealii::hp::QCollection<dim-1>, dealii::hp::QCollection<1>,
             dealii::hp::FECollection<dim> >
-DGBase<dim,real>::create_collection_tuple(const unsigned int max_degree, const int nstate) const
+DGBase<dim,real>::create_collection_tuple(const unsigned int max_degree, const int nstate, const Parameters::AllParameters *const parameters_input) const
 {
     dealii::hp::MappingCollection<dim> mapping_coll;
     dealii::hp::FECollection<dim>      fe_coll;
@@ -137,7 +136,8 @@ DGBase<dim,real>::create_collection_tuple(const unsigned int max_degree, const i
     dealii::hp::QCollection<1>         oned_quad_coll;
 
     dealii::hp::FECollection<dim>      fe_coll_lagr;
-    for (unsigned int degree=0; degree<=max_degree; ++degree) {
+    int minimum_degree = (parameters_input->use_collocated_nodes==true) ?  1 :  0;
+    for (unsigned int degree=minimum_degree; degree<=max_degree; ++degree) {
         //const dealii::MappingQ<dim,dim> mapping(degree, true);
         //const dealii::MappingQ<dim,dim> mapping(degree+1, true);
         const dealii::MappingManifold<dim,dim> mapping;
@@ -147,9 +147,44 @@ DGBase<dim,real>::create_collection_tuple(const unsigned int max_degree, const i
         const dealii::FESystem<dim,dim> fe_system(fe_dg, nstate);
         fe_coll.push_back (fe_system);
 
-        const dealii::QGauss<1>     oned_quad(degree+1);
-        const dealii::QGauss<dim>   volume_quad(degree+1);
-        const dealii::QGauss<dim-1> face_quad(degree+1);
+        //
+
+        dealii::Quadrature<1>     oned_quad(degree+1);
+        dealii::Quadrature<dim>   volume_quad(degree+1);
+        dealii::Quadrature<dim-1> face_quad(degree+1); //removed const
+
+        if (parameters_input->use_collocated_nodes)
+        	{
+        		dealii::QGaussLobatto<1> oned_quad_Gauss_Lobatto (degree+1);
+        		dealii::QGaussLobatto<dim> vol_quad_Gauss_Lobatto (degree+1);
+        		oned_quad = oned_quad_Gauss_Lobatto;
+        		volume_quad = vol_quad_Gauss_Lobatto;
+
+        		if(dim == 1)
+        		{
+        			dealii::QGauss<dim-1> face_quad_Gauss_Legendre (degree+1);
+        			face_quad = face_quad_Gauss_Legendre;
+        		}
+        		else
+        		{
+        			dealii::QGaussLobatto<dim-1> face_quad_Gauss_Lobatto (degree+1);
+        			face_quad = face_quad_Gauss_Lobatto;
+        		}
+
+
+        	}
+        	else
+        	{
+        		dealii::QGauss<1> oned_quad_Gauss_Legendre (degree+1);
+        		dealii::QGauss<dim> vol_quad_Gauss_Legendre (degree+1);
+        		dealii::QGauss<dim-1> face_quad_Gauss_Legendre (degree+1);
+        		oned_quad = oned_quad_Gauss_Legendre;
+        		volume_quad = vol_quad_Gauss_Legendre;
+        		face_quad = face_quad_Gauss_Legendre;
+        	}
+        //
+
+
         volume_quad_coll.push_back (volume_quad);
         face_quad_coll.push_back (face_quad);
         oned_quad_coll.push_back (oned_quad);
@@ -193,7 +228,9 @@ void DGBase<dim,real>::set_triangulation(Triangulation *triangulation_input)
     triangulation = triangulation_input;
     dof_handler.initialize(*triangulation, fe_collection);
 
-    set_all_cells_fe_degree ( max_degree );
+    set_all_cells_fe_degree (fe_collection.size()-1); // Always sets it to the maximum degree
+
+    //set_all_cells_fe_degree ( max_degree-min_degree);
 }
 
 
@@ -252,26 +289,183 @@ void DGBase<dim,real>::assemble_residual (const bool compute_dRdW)
 
             // Case 1:
             // Face at boundary
-            if (current_face->at_boundary()) {
+            if (current_face->at_boundary() && !current_cell->has_periodic_neighbor(iface) ) {
 
                 n_face_visited++;
 
                 fe_values_collection_face_int.reinit (current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell);
-                const dealii::FEFaceValues<dim,dim> &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values();
-                const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1);
-                const unsigned int normal_direction = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
-                const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction);
 
-                real penalty = deg1sq / vol_div_facearea1;
-
-                const unsigned int boundary_id = current_face->boundary_id();
-                // Need to somehow get boundary type from the mesh
-                if ( compute_dRdW ) {
-                    assemble_boundary_term_implicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs);
-                } else { 
-                    assemble_boundary_term_explicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs);
+                if(current_face->at_boundary() && all_parameters->use_periodic_bc == true && dim == 1) //using periodic BCs (for 1d)
+				{
+                	int cell_index  = current_cell->index();
+					//int cell_index = current_cell->index();
+					if (cell_index == 0 && iface == 0)
+					{
+						fe_values_collection_face_int.reinit(current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell);
+						neighbor_cell = dof_handler.begin_active();
+						for (unsigned int i = 0 ; i < triangulation->n_active_cells() - 1; ++i)
+						{
+							++neighbor_cell;
+						}
+						neighbor_cell->get_dof_indices(neighbor_dofs_indices);
+	                 	const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index();
+
+						fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1,fe_index_neigh_cell,fe_index_neigh_cell,fe_index_neigh_cell);
+
+					}
+					else if (cell_index == (int) triangulation->n_active_cells() - 1 && iface == 1)
+					{
+						fe_values_collection_face_int.reinit(current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell);
+						neighbor_cell = dof_handler.begin_active();
+						neighbor_cell->get_dof_indices(neighbor_dofs_indices);
+	                 	const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index();
+						fe_values_collection_face_ext.reinit(neighbor_cell,(iface == 1) ? 0 : 1, fe_index_neigh_cell, fe_index_neigh_cell, fe_index_neigh_cell); //not sure how changing the face number would work in dim!=1-dimensions.
+					}
+
+					//std::cout << "cell " << current_cell->index() << "'s " << iface << "th face has neighbour: " << neighbor_cell->index() << std::endl;
+                	const int neighbor_face_no = (iface ==1) ? 0:1;
+                	const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index();
+
+
+                	const dealii::FEFaceValues<dim,dim> &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values();
+                	const dealii::FEFaceValues<dim,dim> &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values();
+
+                	const dealii::FESystem<dim,dim> &neigh_fe_ref = fe_collection[fe_index_neigh_cell];
+                	const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree();
+                	const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell();
+
+					dealii::Vector<double> neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization
+
+
+                	const unsigned int normal_direction1 = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
+					const unsigned int normal_direction2 = dealii::GeometryInfo<dim>::unit_normal_direction[neighbor_face_no];
+					const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1);
+					const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1);
+
+					//const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children();
+					const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1);
+					const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2);
+
+					const real penalty1 = deg1sq / vol_div_facearea1;
+					const real penalty2 = deg2sq / vol_div_facearea2;
+
+					real penalty = 0.5 * ( penalty1 + penalty2 );
+
+					if ( compute_dRdW ) {
+						assemble_face_term_implicit (
+													fe_values_face_int, fe_values_face_ext,
+													penalty,
+													current_dofs_indices, neighbor_dofs_indices,
+													current_cell_rhs, neighbor_cell_rhs);
+					} else {
+						assemble_face_term_explicit (
+													fe_values_face_int, fe_values_face_ext,
+													penalty,
+													current_dofs_indices, neighbor_dofs_indices,
+													current_cell_rhs, neighbor_cell_rhs);
+					}
+
+				}
+
+                else
+                {
+                	const dealii::FEFaceValues<dim,dim> &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values();
+					const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1);
+					const unsigned int normal_direction = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
+					const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction);
+
+					real penalty = deg1sq / vol_div_facearea1;
+
+					const unsigned int boundary_id = current_face->boundary_id();
+					// Need to somehow get boundary type from the mesh
+					if ( compute_dRdW ) {
+						assemble_boundary_term_implicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs);
+					} else {
+						assemble_boundary_term_explicit (boundary_id, fe_values_face_int, penalty, current_dofs_indices, current_cell_rhs);
+					}
                 }
 
+				//CASE 1.5: periodic boundary conditions
+				//note that periodicity is not adapted for hp adaptivity yet. this needs to be figured out in the future
+			} else if (current_face->at_boundary() && current_cell->has_periodic_neighbor(iface)){
+
+				neighbor_cell = current_cell->periodic_neighbor(iface);
+				//std::cout << "cell " << current_cell->index() << " at boundary" <<std::endl;
+				//std::cout << "periodic neighbour on face " << iface << " is " << neighbor_cell->index() << std::endl;
+
+
+				if (!current_cell->periodic_neighbor_is_coarser(iface) &&
+					(neighbor_cell->index() > current_cell->index() ||
+					 (neighbor_cell->index() == current_cell->index() && current_cell->level() < neighbor_cell->level())
+					)
+				   )
+				{
+					 n_face_visited++;
+					Assert (current_cell->periodic_neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError());
+
+
+					// Corresponding face of the neighbor.
+					// e.g. The 4th face of the current cell might correspond to the 3rd face of the neighbor
+					const unsigned int neighbor_face_no = current_cell->periodic_neighbor_of_periodic_neighbor(iface);
+
+					// Get information about neighbor cell
+					const unsigned int fe_index_neigh_cell = neighbor_cell->active_fe_index();
+					const dealii::FESystem<dim,dim> &neigh_fe_ref = fe_collection[fe_index_neigh_cell];
+					const unsigned int neigh_cell_degree = neigh_fe_ref.tensor_degree();
+					const unsigned int n_dofs_neigh_cell = neigh_fe_ref.n_dofs_per_cell();
+
+					// Local rhs contribution from neighbor
+					dealii::Vector<double> neighbor_cell_rhs (n_dofs_neigh_cell); // Defaults to 0.0 initialization
+
+					// Obtain the mapping from local dof indices to global dof indices for neighbor cell
+					neighbor_dofs_indices.resize(n_dofs_neigh_cell);
+					neighbor_cell->get_dof_indices (neighbor_dofs_indices);
+
+					fe_values_collection_face_int.reinit (current_cell, iface, fe_index_curr_cell, fe_index_curr_cell, fe_index_curr_cell);
+					const dealii::FEFaceValues<dim,dim> &fe_values_face_int = fe_values_collection_face_int.get_present_fe_values();
+					fe_values_collection_face_ext.reinit (neighbor_cell, neighbor_face_no, fe_index_neigh_cell, fe_index_neigh_cell, fe_index_neigh_cell);
+					const dealii::FEFaceValues<dim,dim> &fe_values_face_ext = fe_values_collection_face_ext.get_present_fe_values();
+
+					const unsigned int normal_direction1 = dealii::GeometryInfo<dim>::unit_normal_direction[iface];
+					const unsigned int normal_direction2 = dealii::GeometryInfo<dim>::unit_normal_direction[neighbor_face_no];
+					const unsigned int deg1sq = (curr_cell_degree == 0) ? 1 : curr_cell_degree * (curr_cell_degree+1);
+					const unsigned int deg2sq = (neigh_cell_degree == 0) ? 1 : neigh_cell_degree * (neigh_cell_degree+1);
+
+					//const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1) / current_face->number_of_children();
+					const real vol_div_facearea1 = current_cell->extent_in_direction(normal_direction1);
+					const real vol_div_facearea2 = neighbor_cell->extent_in_direction(normal_direction2);
+
+					const real penalty1 = deg1sq / vol_div_facearea1;
+					const real penalty2 = deg2sq / vol_div_facearea2;
+
+					real penalty = 0.5 * ( penalty1 + penalty2 );
+					//penalty = 1;//99;
+
+					if ( compute_dRdW ) {
+						assemble_face_term_implicit (
+								fe_values_face_int, fe_values_face_ext,
+								penalty,
+								current_dofs_indices, neighbor_dofs_indices,
+								current_cell_rhs, neighbor_cell_rhs);
+					} else {
+						assemble_face_term_explicit (
+								fe_values_face_int, fe_values_face_ext,
+								penalty,
+								current_dofs_indices, neighbor_dofs_indices,
+								current_cell_rhs, neighbor_cell_rhs);
+					}
+
+					// Add local contribution from neighbor cell to global vector
+					for (unsigned int i=0; i<n_dofs_neigh_cell; ++i) {
+						right_hand_side(neighbor_dofs_indices[i]) += neighbor_cell_rhs(i);
+					}
+				}
+				else
+				{
+					//do nothing
+				}
+
+
             // Case 2:
             // Neighbour is finer occurs if the face has children
             // In this case, we loop over the current large face's subfaces and visit multiple neighbors
@@ -335,7 +529,6 @@ void DGBase<dim,real>::assemble_residual (const bool compute_dRdW)
                             current_dofs_indices, neighbor_dofs_indices,
                             current_cell_rhs, neighbor_cell_rhs);
                     }
-
                     // Add local contribution from neighbor cell to global vector
                     for (unsigned int i=0; i<n_dofs_neigh_cell; ++i) {
                         right_hand_side(neighbor_dofs_indices[i]) += neighbor_cell_rhs(i);
@@ -345,8 +538,8 @@ void DGBase<dim,real>::assemble_residual (const bool compute_dRdW)
             // Case 3:
             // Neighbor cell is NOT coarser
             // Therefore, they have the same coarseness, and we need to choose one of them to do the work
-            } else if (
-                !current_cell->neighbor_is_coarser(iface) &&
+            } else if ( //added a criteria for periodicity
+                (!current_cell->neighbor_is_coarser(iface))  &&
                     // Cell with lower index does work
                     (neighbor_cell->index() > current_cell->index() || 
                     // If both cells have same index
@@ -358,7 +551,7 @@ void DGBase<dim,real>::assemble_residual (const bool compute_dRdW)
                 n_face_visited++;
                 Assert (current_cell->neighbor(iface).state() == dealii::IteratorState::valid, dealii::ExcInternalError());
 
-                auto neighbor_cell = current_cell->neighbor(iface);
+                auto neighbor_cell = current_cell->neighbor_or_periodic_neighbor(iface);
                 // Corresponding face of the neighbor.
                 // e.g. The 4th face of the current cell might correspond to the 3rd face of the neighbor
                 const unsigned int neighbor_face_no = current_cell->neighbor_of_neighbor(iface);
@@ -447,6 +640,7 @@ unsigned int DGBase<dim,real>::n_dofs () const
 template <int dim, typename real>
 void DGBase<dim,real>::output_results_vtk (const unsigned int cycle)// const
 {
+
     dealii::DataOut<dim, dealii::hp::DoFHandler<dim>> data_out;
     data_out.attach_dof_handler (dof_handler);
 
@@ -471,11 +665,19 @@ void DGBase<dim,real>::output_results_vtk (const unsigned int cycle)// const
     std::vector<unsigned int> active_fe_indices;
     dof_handler.get_active_fe_indices(active_fe_indices);
     dealii::Vector<double> active_fe_indices_dealiivector(active_fe_indices.begin(), active_fe_indices.end());
-    //using DVTenum = dealii::DataOut_DoFData<dealii::hp::DoFHandler<dim>,dim>::DataVectorType;
+    dealii::Vector<double> cell_poly_degree = active_fe_indices_dealiivector;
+
+//    int index = 0;
+//    for (auto current_cell_poly = cell_poly_degree.begin(); current_cell_poly != cell_poly_degree.end(); ++current_cell_poly) {
+//    	current_cell_poly[index] = fe_collection[active_fe_indices_dealiivector[index]].tensor_degree();
+//    	index++;
+//    }
+//    //using DVTenum = dealii::DataOut_DoFData<dealii::hp::DoFHandler<dim>,dim>::DataVectorType;
+//    data_out.add_data_vector (cell_poly_degree, "PolynomialDegree", dealii::DataOut_DoFData<dealii::hp::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
     data_out.add_data_vector (active_fe_indices_dealiivector, "PolynomialDegree", dealii::DataOut_DoFData<dealii::hp::DoFHandler<dim>,dim>::DataVectorType::type_cell_data);
 
 
-    assemble_residual (false);
+    //assemble_residual (false);
     std::vector<std::string> residual_names;
     for(int s=0;s<nstate;++s) {
         std::string varname = "residual" + dealii::Utilities::int_to_string(s,1);
diff --git a/src/dg/dg.h b/src/dg/dg.h
index 11f35fe4e..6dfed5e17 100644
--- a/src/dg/dg.h
+++ b/src/dg/dg.h
@@ -61,6 +61,8 @@ class DGBase
      *  DGBase cannot use nstate as a compile-time known.  */
     const unsigned int max_degree;
 
+
+
     /// Pointer to all parameters
     const Parameters::AllParameters *const all_parameters;
 
@@ -249,7 +251,7 @@ class DGBase
     std::tuple< dealii::hp::MappingCollection<dim>, dealii::hp::FECollection<dim>,
                 dealii::hp::QCollection<dim>, dealii::hp::QCollection<dim-1>, dealii::hp::QCollection<1>,
                 dealii::hp::FECollection<dim> >
-                create_collection_tuple(const unsigned int max_degree, const int nstate) const;
+                create_collection_tuple(const unsigned int max_degree, const int nstate, const Parameters::AllParameters *const parameters_input) const;
 
     /// Evaluate the integral over the cell volume
     virtual void assemble_volume_terms_implicit(
@@ -368,6 +370,8 @@ class DGWeak : public DGBase<dim, real>
     ~DGWeak();
 
 private:
+
+
     /// Contains the physics of the PDE
     std::shared_ptr < Physics::PhysicsBase<dim, nstate, Sacado::Fad::DFad<real> > > pde_physics;
     /// Convective numerical flux
@@ -454,6 +458,7 @@ class DGStrong : public DGBase<dim, real>
     /// Dissipative numerical flux
     NumericalFlux::NumericalFluxDissipative<dim, nstate, Sacado::Fad::DFad<real> > *diss_num_flux;
 
+
     /// Contains the physics of the PDE
     std::shared_ptr < Physics::PhysicsBase<dim, nstate, real > > pde_physics_double;
     /// Convective numerical flux
@@ -461,7 +466,6 @@ class DGStrong : public DGBase<dim, real>
     /// Dissipative numerical flux
     NumericalFlux::NumericalFluxDissipative<dim, nstate, real > *diss_num_flux_double;
 
-
     /// Evaluate the integral over the cell volume
     void assemble_volume_terms_implicit(
         const dealii::FEValues<dim,dim> &fe_values_volume,
diff --git a/src/dg/strong_dg.cpp b/src/dg/strong_dg.cpp
index 328a62720..76aca2e12 100644
--- a/src/dg/strong_dg.cpp
+++ b/src/dg/strong_dg.cpp
@@ -50,6 +50,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
     const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
     dealii::Vector<real> &local_rhs_int_cell)
 {
+
     using ADtype = Sacado::Fad::DFad<real>;
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
@@ -98,6 +99,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
         // Evaluate physical convective flux and source term
         conv_phys_flux_at_q[iquad] = pde_physics->convective_flux (soln_at_q[iquad]);
         diss_phys_flux_at_q[iquad] = pde_physics->dissipative_flux (soln_at_q[iquad], soln_grad_at_q[iquad]);
+
         if(this->all_parameters->manufactured_convergence_study_param.use_manufactured_source_term) {
             source_at_q[iquad] = pde_physics->source_term (fe_values_vol.quadrature_point(iquad), soln_at_q[iquad]);
         }
@@ -108,12 +110,17 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
     // Since we have nodal values of the flux, we use the Lagrange polynomials to obtain the gradients at the quadrature points.
     const dealii::FEValues<dim,dim> &fe_values_lagrange = this->fe_values_collection_volume_lagrange.get_present_fe_values();
     std::vector<ADArray> flux_divergence(n_quad_pts);
+
+    std::array<std::array<std::vector<ADtype>,nstate>,dim> f;
+    std::array<std::array<std::vector<ADtype>,nstate>,dim> g;
+
     for (int istate = 0; istate<nstate; ++istate) {
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
             flux_divergence[iquad][istate] = 0.0;
             for ( unsigned int flux_basis = 0; flux_basis < n_quad_pts; ++flux_basis ) {
                 flux_divergence[iquad][istate] += conv_phys_flux_at_q[flux_basis][istate] * fe_values_lagrange.shape_grad(flux_basis,iquad);
             }
+
         }
     }
 
@@ -129,6 +136,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
 
         ADtype rhs = 0;
 
+
         const unsigned int istate = fe_values_vol.get_fe().system_to_component_index(itest).first;
 
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
@@ -136,6 +144,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_implicit(
             // Convective
             // Now minus such 2 integrations by parts
             assert(JxW[iquad] - fe_values_lagrange.JxW(iquad) < 1e-14);
+
             rhs = rhs - fe_values_vol.shape_value_component(itest,iquad,istate) * flux_divergence[iquad][istate] * JxW[iquad];
 
             //// Diffusive
@@ -470,6 +479,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
     const std::vector<dealii::types::global_dof_index> &cell_dofs_indices,
     dealii::Vector<real> &local_rhs_int_cell)
 {
+	//std::cout << "assembling cell terms" << std::endl;
     using ADtype = Sacado::Fad::DFad<real>;
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
@@ -532,7 +542,14 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
             flux_divergence[iquad][istate] = 0.0;
             for ( unsigned int flux_basis = 0; flux_basis < n_quad_pts; ++flux_basis ) {
-                flux_divergence[iquad][istate] += conv_phys_flux_at_q[flux_basis][istate] * fe_values_lagrange.shape_grad(flux_basis,iquad);
+            	if (this->all_parameters->use_split_form == true)
+            	{
+            		flux_divergence[iquad][istate] += 2* pde_physics->convective_numerical_split_flux(soln_at_q[iquad],soln_at_q[flux_basis])[istate] *  fe_values_lagrange.shape_grad(flux_basis,iquad);
+            	}
+            	else
+            	{
+                    flux_divergence[iquad][istate] += conv_phys_flux_at_q[flux_basis][istate] * fe_values_lagrange.shape_grad(flux_basis,iquad);
+            	}
             }
         }
     }
@@ -556,6 +573,7 @@ void DGStrong<dim,nstate,real>::assemble_volume_terms_explicit(
             // Convective
             // Now minus such 2 integrations by parts
             assert(JxW[iquad] - fe_values_lagrange.JxW(iquad) < 1e-14);
+
             rhs = rhs - fe_values_vol.shape_value_component(itest,iquad,istate) * flux_divergence[iquad][istate] * JxW[iquad];
 
             //// Diffusive
@@ -720,6 +738,7 @@ void DGStrong<dim,nstate,real>::assemble_face_term_explicit(
     dealii::Vector<real>          &local_rhs_int_cell,
     dealii::Vector<real>          &local_rhs_ext_cell)
 {
+	//std::cout << "assembling face terms" << std::endl;
     using ADtype = Sacado::Fad::DFad<real>;
     using ADArray = std::array<ADtype,nstate>;
     using ADArrayTensor1 = std::array< dealii::Tensor<1,dim,ADtype>, nstate >;
@@ -808,11 +827,16 @@ void DGStrong<dim,nstate,real>::assemble_face_term_explicit(
         //std::cout << "Energy ext" << soln_ext[iquad][nstate-1] << std::endl;
 
         // Evaluate physical convective flux, physical dissipative flux, and source term
+
+        //std::cout <<"evaluating numerical fluxes" <<std::endl;
         conv_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
 
         conv_phys_flux_int[iquad] = pde_physics->convective_flux (soln_int[iquad]);
         conv_phys_flux_ext[iquad] = pde_physics->convective_flux (soln_ext[iquad]);
 
+       // std::cout <<"done evaluating numerical fluxes" <<std::endl;
+
+
         diss_soln_num_flux[iquad] = diss_num_flux->evaluate_solution_flux(soln_int[iquad], soln_ext[iquad], normal_int);
 
         ADArrayTensor1 diss_soln_jump_int, diss_soln_jump_ext;
diff --git a/src/numerical_flux/CMakeLists.txt b/src/numerical_flux/CMakeLists.txt
index 67e387273..f105464d2 100644
--- a/src/numerical_flux/CMakeLists.txt
+++ b/src/numerical_flux/CMakeLists.txt
@@ -1,6 +1,7 @@
 set(FLUX_SOURCE
     viscous_numerical_flux.cpp
     numerical_flux.cpp
+    split_form_numerical_flux.cpp
     )
 
 foreach(dim RANGE 1 3)
diff --git a/src/numerical_flux/numerical_flux.cpp b/src/numerical_flux/numerical_flux.cpp
index 115ef2dae..999def696 100644
--- a/src/numerical_flux/numerical_flux.cpp
+++ b/src/numerical_flux/numerical_flux.cpp
@@ -1,7 +1,8 @@
 #include <Sacado.hpp>
 #include <deal.II/differentiation/ad/sacado_product_types.h>
 #include "numerical_flux.h"
-
+#include "viscous_numerical_flux.h"
+#include "split_form_numerical_flux.h"
 
 namespace PHiLiP {
 namespace NumericalFlux {
@@ -34,6 +35,9 @@ ::create_convective_numerical_flux(
     } else if(conv_num_flux_type == AllParam::roe) {
         if constexpr (dim+2==nstate) return new Roe<dim, nstate, real>(physics_input);
     }
+    else if (conv_num_flux_type == AllParam::split_form) {
+    	return new SplitFormNumFlux<dim, nstate, real>(physics_input);
+    }
 
     std::cout << "Invalid numerical flux" << std::endl;
     return nullptr;
diff --git a/src/numerical_flux/split_form_numerical_flux.cpp b/src/numerical_flux/split_form_numerical_flux.cpp
new file mode 100644
index 000000000..dfa30347b
--- /dev/null
+++ b/src/numerical_flux/split_form_numerical_flux.cpp
@@ -0,0 +1,49 @@
+#include "split_form_numerical_flux.h"
+#include <Sacado.hpp>
+#include <deal.II/differentiation/ad/sacado_product_types.h>
+
+namespace PHiLiP {
+namespace NumericalFlux {
+
+template <int dim, int nstate, typename real>
+std::array<real, nstate> SplitFormNumFlux<dim,nstate,real>::evaluate_flux(
+	const std::array<real, nstate> &soln_int,
+    const std::array<real, nstate> &soln_ext,
+    const dealii::Tensor<1,dim,real> &normal_int) const
+	{
+		//std::cout << "evaluating the split form flux" <<std::endl;
+		using RealArrayVector = std::array<dealii::Tensor<1,dim,real>,nstate>;
+	    RealArrayVector conv_phys_split_flux;
+
+	    conv_phys_split_flux = pde_physics->convective_numerical_split_flux (soln_int,soln_ext);
+	    //std::cout << "done evaluating the conv num split flux" <<std::endl;
+
+	    const real conv_max_eig_int = pde_physics->max_convective_eigenvalue(soln_int);
+	   // std::cout << "1st eig" << std::endl;
+	    const real conv_max_eig_ext = pde_physics->max_convective_eigenvalue(soln_ext);
+	    //std::cout << "2nd eig" << std::endl;
+	    const real conv_max_eig = std::max(conv_max_eig_int, conv_max_eig_ext);
+
+	   // std::cout << "obtained the max eig" <<std::endl;
+	    // Scalar dissipation
+	    std::array<real, nstate> numerical_flux_dot_n;
+	    for (int s=0; s<nstate; s++) {
+	        numerical_flux_dot_n[s] = conv_phys_split_flux[s]*normal_int - 0.5 * conv_max_eig * (soln_ext[s]-soln_int[s]);
+	    }
+	   // std::cout << "about to return split num flux" <<std::endl;
+	    return numerical_flux_dot_n;
+	}
+
+template class SplitFormNumFlux<PHILIP_DIM, 1, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 2, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 2, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 3, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 3, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 4, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 4, Sacado::Fad::DFad<double> >;
+template class SplitFormNumFlux<PHILIP_DIM, 5, double>;
+template class SplitFormNumFlux<PHILIP_DIM, 5, Sacado::Fad::DFad<double> >;
+
+}
+}
diff --git a/src/numerical_flux/split_form_numerical_flux.h b/src/numerical_flux/split_form_numerical_flux.h
new file mode 100644
index 000000000..22116c5dd
--- /dev/null
+++ b/src/numerical_flux/split_form_numerical_flux.h
@@ -0,0 +1,40 @@
+#ifndef __SPLIT_NUM_FLUX__
+#define __SPLIT_NUM_FLUX__
+
+#include <deal.II/base/tensor.h>
+#include "physics/physics.h"
+#include "numerical_flux/numerical_flux.h"
+
+namespace PHiLiP {
+namespace NumericalFlux {
+
+/// Lax-Friedrichs numerical flux. Derived from NumericalFluxConvective.
+template<int dim, int nstate, typename real>
+class SplitFormNumFlux: public NumericalFluxConvective<dim, nstate, real>
+{
+public:
+
+/// Constructor
+SplitFormNumFlux(std::shared_ptr <Physics::PhysicsBase<dim, nstate, real>> physics_input)
+:
+pde_physics(physics_input)
+{};
+/// Destructor
+~SplitFormNumFlux() {};
+
+/// Returns the convective numerical flux at an interface.
+std::array<real, nstate> evaluate_flux (
+    const std::array<real, nstate> &soln_int,
+    const std::array<real, nstate> &soln_ext,
+    const dealii::Tensor<1,dim,real> &normal1) const;
+
+protected:
+/// Numerical flux requires physics to evaluate convective eigenvalues.
+const std::shared_ptr < Physics::PhysicsBase<dim, nstate, real> > pde_physics;
+
+};
+
+}
+}
+
+#endif
diff --git a/src/ode_solver/ode_solver.cpp b/src/ode_solver/ode_solver.cpp
index 492a38eab..961a7e439 100644
--- a/src/ode_solver/ode_solver.cpp
+++ b/src/ode_solver/ode_solver.cpp
@@ -96,8 +96,9 @@ int ODESolver<dim,real>::steady_state ()
                   << std::endl;
 
         if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-            (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-        pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+            (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+		pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+	}
 
         this->dg->assemble_residual ();
         this->residual_norm = this->dg->get_residual_l2norm();
@@ -120,7 +121,6 @@ int ODESolver<dim,real>::steady_state ()
                 this->dg->output_results_vtk(file_number);
             }
         }
-
     }
     return 1;
 }
@@ -146,21 +146,29 @@ int ODESolver<dim,real>::advance_solution_time (double time_advance)
     while (this->current_iteration < number_of_time_steps)
     {
         if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-            (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-        pcout << " ********************************************************** "
-                  << std::endl
-                  << " Iteration: " << this->current_iteration + 1
-                  << " out of: " << number_of_time_steps
-                  << std::endl;
+            (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+		pcout << " ********************************************************** "
+			  << std::endl
+			  << " Iteration: " << this->current_iteration + 1
+			  << " out of: " << number_of_time_steps
+			  << std::endl;
+	}
+        dg->assemble_residual(false);
 
         if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-            (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-        pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+            (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+		pcout << " Evaluating right-hand side and setting system_matrix to Jacobian... " << std::endl;
+	}
 
         step_in_time(constant_time_step);
 
+
+   if (this->current_iteration%ode_param.print_iteration_modulo == 0) {
+		this->dg->output_results_vtk(this->current_iteration);
+	}
         ++(this->current_iteration);
 
+
         //this->dg->output_results_vtk(this->current_iteration);
     }
     return 1;
@@ -181,8 +189,9 @@ void Implicit_ODESolver<dim,real>::step_in_time (real dt)
     this->dg->add_mass_matrices(1.0/dt);
 
     if ((ode_param.ode_output) == Parameters::OutputEnum::verbose &&
-        (this->current_iteration%ode_param.print_iteration_modulo) == 0 )
-    pcout << " Evaluating system update... " << std::endl;
+        (this->current_iteration%ode_param.print_iteration_modulo) == 0 ) {
+	    pcout << " Evaluating system update... " << std::endl;
+    }
 
     solve_linear (
         this->dg->system_matrix,
@@ -200,7 +209,7 @@ void Explicit_ODESolver<dim,real>::step_in_time (real dt)
 {
     // this->dg->assemble_residual (); // Not needed since it is called in the base class for time step
     this->current_time += dt;
-    const int rk_order = 3;
+    const int rk_order = 1;
     if (rk_order == 1) {
         this->dg->global_inverse_mass_matrix.vmult(this->solution_update, this->dg->right_hand_side);
         this->update_norm = this->solution_update.l2_norm();
diff --git a/src/parameters/all_parameters.cpp b/src/parameters/all_parameters.cpp
index b4bf4a163..f3181d96d 100644
--- a/src/parameters/all_parameters.cpp
+++ b/src/parameters/all_parameters.cpp
@@ -25,6 +25,19 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm)
     prm.declare_entry("use_weak_form", "true",
                       dealii::Patterns::Bool(),
                       "Use weak form by default. If false, use strong form.");
+
+    prm.declare_entry("use_collocated_nodes", "false",
+    				  dealii::Patterns::Bool(),
+					  "Use Gauss-Legendre by default. Otherwise, use Gauss-Lobatto to collocate.");
+
+    prm.declare_entry("use_split_form", "false",
+    				  dealii::Patterns::Bool(),
+					  "Use original form by defualt. Otherwise, split the fluxes.");
+
+    prm.declare_entry("use_periodic_bc", "false",
+    				  dealii::Patterns::Bool(),
+					  "Use other boundary conditions by default. Otherwise use periodic (for 1d burgers only");
+
     prm.declare_entry("test_type", "run_control",
                       dealii::Patterns::Selection(
                       " run_control | "
@@ -51,7 +64,8 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm)
                           " convection_diffusion | "
                           " advection_vector | "
                           " burgers_inviscid | "
-                          " euler"),
+                          " euler |"
+                    	  " mhd"),
                       "The PDE we want to solve. "
                       "Choices are " 
                       " <advection | " 
@@ -59,11 +73,13 @@ void AllParameters::declare_parameters (dealii::ParameterHandler &prm)
                       "  convection_diffusion | "
                       "  advection_vector | "
                       "  burgers_inviscid | "
-                      "  euler>.");
+                      "  euler | "
+				      "  mhd>.");
     prm.declare_entry("conv_num_flux", "lax_friedrichs",
-                      dealii::Patterns::Selection("lax_friedrichs | roe"),
+                      dealii::Patterns::Selection("lax_friedrichs | roe | split_form"),
                       "Convective numerical flux. "
-                      "Choices are <lax_friedrichs | roe>.");
+                      "Choices are <lax_friedrichs | roe | split_form>.");
+
     prm.declare_entry("diss_num_flux", "symm_internal_penalty",
                       dealii::Patterns::Selection("symm_internal_penalty"),
                       "Dissipative numerical flux. "
@@ -113,10 +129,15 @@ void AllParameters::parse_parameters (dealii::ParameterHandler &prm)
         pde_type = euler;
         nstate = dimension+2;
     }
+
     use_weak_form = prm.get_bool("use_weak_form");
+    use_collocated_nodes = prm.get_bool("use_collocated_nodes");
+    use_split_form = prm.get_bool("use_split_form");
+    use_periodic_bc = prm.get_bool("use_periodic_bc");
 
     const std::string conv_num_flux_string = prm.get("conv_num_flux");
     if (conv_num_flux_string == "lax_friedrichs") conv_num_flux_type = lax_friedrichs;
+    if (conv_num_flux_string == "split_form") conv_num_flux_type = split_form;
     if (conv_num_flux_string == "roe") conv_num_flux_type = roe;
 
     const std::string diss_num_flux_string = prm.get("diss_num_flux");
diff --git a/src/parameters/all_parameters.h b/src/parameters/all_parameters.h
index 019ad2fe0..606344976 100644
--- a/src/parameters/all_parameters.h
+++ b/src/parameters/all_parameters.h
@@ -35,6 +35,14 @@ class AllParameters
     /// Flag to use weak or strong form of DG
     bool use_weak_form;
 
+    /// Flag to use Gauss-Lobatto Nodes;
+    bool use_collocated_nodes;
+
+    /// Flag to use split form.
+    bool use_split_form;
+
+    bool use_periodic_bc;
+
     /// Number of state variables. Will depend on PDE
     int nstate;
 
@@ -45,8 +53,12 @@ class AllParameters
         euler_cylinder,
         euler_vortex,
         euler_entropy_waves,
+		euler_split_taylor_green,
         numerical_flux_convervation,
-        jacobian_regression};
+        jacobian_regression,
+    	burgers_split_form,
+    	advection_periodicity,
+    	};
     TestType test_type;
 
     /// Currently allows to solve advection, diffusion, convection-diffusion
@@ -56,7 +68,8 @@ class AllParameters
         convection_diffusion,
         advection_vector,
         burgers_inviscid,
-        euler};
+        euler,
+		mhd};
 
     /// Possible boundary types, NOT IMPLEMENTED YET
     enum BoundaryType { 
@@ -74,8 +87,10 @@ class AllParameters
     /// Store the PDE type to be solved
     PartialDifferentialEquation pde_type;
 
-    /// Currently only Lax-Friedrichs and roe can be used as an input parameter
-    enum ConvectiveNumericalFlux { lax_friedrichs, roe };
+
+    /// Currently only Lax-Friedrichs, roe, and split_form can be used as an input parameter
+    enum ConvectiveNumericalFlux { lax_friedrichs, roe, split_form};
+
     /// Store convective flux type
     ConvectiveNumericalFlux conv_num_flux_type;
 
diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt
index aae1c5376..91f987ff4 100644
--- a/src/physics/CMakeLists.txt
+++ b/src/physics/CMakeLists.txt
@@ -5,6 +5,7 @@ SET(SOURCE
     burgers.cpp
     euler.cpp
     manufactured_solution.cpp
+    # mhd.cpp
     )
 
 foreach(dim RANGE 1 3)
diff --git a/src/physics/burgers.cpp b/src/physics/burgers.cpp
index e98308292..7b9d97ac4 100644
--- a/src/physics/burgers.cpp
+++ b/src/physics/burgers.cpp
@@ -69,6 +69,20 @@ ::convective_flux (const std::array<real,nstate> &solution) const
     return conv_flux;
 }
 
+template <int dim, int nstate, typename real>
+std::array<dealii::Tensor<1,dim,real>,nstate> Burgers<dim,nstate,real>::convective_numerical_split_flux (
+        		const std::array<real,nstate> &soln_const,
+				const std::array<real,nstate> & soln_loop) const
+{
+	std::array<dealii::Tensor<1,dim,real>,nstate> conv_flux;
+	    for (int flux_dim=0; flux_dim<dim; ++flux_dim) {
+	        for (int s=0; s<nstate; ++s) {
+	            conv_flux[s][flux_dim] = 1./6. * (soln_const[flux_dim]*soln_const[flux_dim] + soln_const[flux_dim]*soln_loop[s] + soln_loop[s]*soln_loop[s]);
+	        }
+	    }
+	    return conv_flux;
+}
+
 template <int dim, int nstate, typename real>
 real Burgers<dim,nstate,real>
 ::diffusion_coefficient () const
diff --git a/src/physics/burgers.h b/src/physics/burgers.h
index 14b725168..7c120e156 100644
--- a/src/physics/burgers.h
+++ b/src/physics/burgers.h
@@ -43,6 +43,11 @@ class Burgers : public PhysicsBase <dim, nstate, real>
     /// Convective flux: \f$ \mathbf{F}_{conv} =  u \f$
     std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (const std::array<real,nstate> &solution) const;
 
+    /// Convective split flux
+    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux (
+        		const std::array<real,nstate> &soln_const,
+				const std::array<real,nstate> & soln_loop) const;
+
     /// Spectral radius of convective term Jacobian is 'c'
     std::array<real,nstate> convective_eigenvalues (
         const std::array<real,nstate> &/*solution*/,
diff --git a/src/physics/convection_diffusion.h b/src/physics/convection_diffusion.h
index ec3e5df45..dea5b2483 100644
--- a/src/physics/convection_diffusion.h
+++ b/src/physics/convection_diffusion.h
@@ -42,6 +42,18 @@ class ConvectionDiffusion : public PhysicsBase <dim, nstate, real>
     /// Convective flux: \f$ \mathbf{F}_{conv} =  u \f$
     std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (const std::array<real,nstate> &solution) const;
 
+
+    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux (
+       		const std::array<real,nstate> &soln_const, const std::array<real,nstate> &soln_loop) const
+		{
+    	std::array<real,nstate> arr_avg;
+    		for (int i = 0 ; i < nstate; ++i)
+    		{
+    			arr_avg[i] = (soln_const[i] + soln_loop[i])/2.;
+    		}
+    		return convective_flux(arr_avg);
+		};
+
     /// Spectral radius of convective term Jacobian is 'c'
     std::array<real,nstate> convective_eigenvalues (
         const std::array<real,nstate> &/*solution*/,
diff --git a/src/physics/euler.cpp b/src/physics/euler.cpp
index 10f9ef56d..8d6a9b121 100644
--- a/src/physics/euler.cpp
+++ b/src/physics/euler.cpp
@@ -1,4 +1,4 @@
-#include <cmath>
+ #include <cmath>
 #include <vector>
 
 #include <Sacado.hpp>
@@ -192,10 +192,20 @@ inline real Euler<dim,nstate,real>
 ::compute_pressure ( const std::array<real,nstate> &conservative_soln ) const
 {
     const real density = conservative_soln[0];
+    //std::cout << "density " << density << std::endl;
+
     const real tot_energy  = conservative_soln[nstate-1];
+  //  std::cout << "tot_energy " << tot_energy << std::endl;
+
     const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+   // std::cout << "vel1 " << vel[0] << std::endl
+   // 					 << vel[1] << std::endl
+//						 << vel[2] <<std::endl;
+
     const real vel2 = compute_velocity_squared(vel);
+    //std::cout << "vel ^2 " << vel2 <<std::endl;
     real pressure = gamm1*(tot_energy - 0.5*density*vel2);
+    //std::cout << "calculated pressure is" << pressure << std::endl;
     if(pressure<0.0) {
         std::cout<<"Cannot compute pressure..."<<std::endl;
         std::cout<<"density "<<density<<std::endl;
@@ -219,7 +229,9 @@ ::compute_sound ( const std::array<real,nstate> &conservative_soln ) const
     }
     assert(density > 0);
     const real pressure = compute_pressure(conservative_soln);
+    //std::cout << "pressure is" << pressure << std::endl;
     const real sound = std::sqrt(pressure*gam/density);
+    //std::cout << "sound is " << sound << std::endl;
     return sound;
 }
 
@@ -243,6 +255,45 @@ ::compute_mach_number ( const std::array<real,nstate> &conservative_soln ) const
     return mach_number;
 }
 
+// Split form functions:
+
+template <int dim, int nstate, typename real>
+inline real Euler<dim,nstate,real>::
+compute_mean_density(const std::array<real,nstate> &soln_const,
+						  const std::array<real,nstate> &soln_loop) const
+{
+	return (soln_const[0] + soln_loop[0])/2.;
+}
+
+template <int dim, int nstate, typename real>
+inline real Euler<dim,nstate,real>::
+compute_mean_pressure(const std::array<real,nstate> &soln_const,
+    				  const std::array<real,nstate> &soln_loop) const
+{
+	real pressure_const = compute_pressure(soln_const);
+	real pressure_loop = compute_pressure(soln_loop);
+	return (pressure_const + pressure_loop)/2.;
+}
+
+template <int dim, int nstate, typename real>
+inline dealii::Tensor<1,dim,real> Euler<dim,nstate,real>::
+compute_mean_velocities(const std::array<real,nstate> &soln_const,
+    				    const std::array<real,nstate> &soln_loop) const
+{
+	dealii::Tensor<1,dim,real> vel_const = compute_velocities(soln_const);
+	dealii::Tensor<1,dim,real> vel_loop = compute_velocities(soln_loop);
+	return (vel_const + vel_loop)/2.;
+}
+
+template <int dim, int nstate, typename real>
+inline real Euler<dim,nstate,real>::
+compute_mean_specific_energy(const std::array<real,nstate> &soln_const,
+    						 const std::array<real,nstate> &soln_loop) const
+{
+	return ((soln_const[nstate-1]/soln_const[0]) + (soln_loop[nstate-1]/soln_loop[0]))/2.;
+}
+
+
 template <int dim, int nstate, typename real>
 std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim,nstate,real>
 ::convective_flux (const std::array<real,nstate> &conservative_soln) const
@@ -268,6 +319,33 @@ ::convective_flux (const std::array<real,nstate> &conservative_soln) const
     return conv_flux;
 }
 
+template <int dim, int nstate, typename real>
+std::array<dealii::Tensor<1,dim,real>,nstate> Euler<dim, nstate, real>
+::convective_numerical_split_flux(const std::array<real,nstate> &soln_const,
+								  const std::array<real,nstate> &soln_loop) const
+{
+	std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
+	const real mean_density = compute_mean_density(soln_const, soln_loop);
+	const real mean_pressure = compute_mean_pressure(soln_const, soln_loop);
+	const dealii::Tensor<1,dim,real> mean_velocities = compute_mean_velocities(soln_const,soln_loop);
+	const real mean_specific_energy = compute_mean_specific_energy(soln_const, soln_loop);
+
+	for (int flux_dim = 0; flux_dim < dim; ++flux_dim)
+	{
+		// Density equation
+		conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];//conservative_soln[1+flux_dim];
+		// Momentum equation
+		for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
+			conv_num_split_flux[1+velocity_dim][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_velocities[velocity_dim];
+		}
+		conv_num_split_flux[1+flux_dim][flux_dim] += mean_pressure; // Add diagonal of pressure
+		// Energy equation
+		conv_num_split_flux[nstate-1][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_specific_energy + mean_pressure * mean_velocities[flux_dim];
+	}
+
+	return conv_num_split_flux;
+}
+
 template <int dim, int nstate, typename real>
 std::array<real,nstate> Euler<dim,nstate,real>
 ::convective_normal_flux (const std::array<real,nstate> &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const
@@ -360,10 +438,22 @@ template <int dim, int nstate, typename real>
 real Euler<dim,nstate,real>
 ::max_convective_eigenvalue (const std::array<real,nstate> &conservative_soln) const
 {
+	//std::cout << "going to calculate max eig" << std::endl;
     const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+	//std::cout << "velocities calculated" << std::endl;
+
     const real sound = compute_sound (conservative_soln);
-    const real vel2 = compute_velocity_squared(vel);
+	//std::cout << "sound calculated" << std::endl;
+
+    /*const*/ real vel2 = compute_velocity_squared(vel);
+	//std::cout << "vel2 calculated" << std::endl;
+
+	if (vel2 < 0.0001)
+		vel2 = 0.0001;
+
     const real max_eig = sqrt(vel2) + sound;
+	//std::cout << "max eig calculated" << std::endl;
+
     return max_eig;
 }
 
diff --git a/src/physics/euler.h b/src/physics/euler.h
index a94302497..b71c9a6c7 100644
--- a/src/physics/euler.h
+++ b/src/physics/euler.h
@@ -155,6 +155,10 @@ class Euler : public PhysicsBase <dim, nstate, real>
     std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (
         const std::array<real,nstate> &conservative_soln) const;
 
+    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux (
+           		const std::array<real,nstate> &soln_const, const std::array<real,nstate> &soln_loop) const;
+
+
     std::array<real,nstate> convective_normal_flux (const std::array<real,nstate> &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const;
 
     /// Convective flux Jacobian: \f$ \frac{\partial \mathbf{F}_{conv}}{\partial w} \cdot \mathbf{n} \f$
@@ -242,6 +246,20 @@ class Euler : public PhysicsBase <dim, nstate, real>
     /** See the book I do like CFD, sec 4.14.2 */
     real compute_temperature_from_density_pressure ( const real density, const real pressure ) const;
 
+    ///These functions are only relevant to the split form. The Euler split form is that of Kennedy & Gruber.
+    /** Refer to Gassner's paper for more information:  */
+    real compute_mean_density(const std::array<real,nstate> &soln_const,
+    						  const std::array<real,nstate> &soln_loop) const;
+
+    real compute_mean_pressure(const std::array<real,nstate> &soln_const,
+        					   const std::array<real,nstate> &soln_loop) const;
+
+    dealii::Tensor<1,dim,real> compute_mean_velocities(const std::array<real,nstate> &soln_const,
+        						  	  	  	  	  	   const std::array<real,nstate> &soln_loop) const;
+
+    real compute_mean_specific_energy(const std::array<real,nstate> &soln_const,
+        						  const std::array<real,nstate> &soln_loop) const;
+
     void boundary_face_values (
         const int /*boundary_type*/,
         const dealii::Point<dim, double> &/*pos*/,
diff --git a/src/physics/mhd.cpp b/src/physics/mhd.cpp
new file mode 100644
index 000000000..df843ee03
--- /dev/null
+++ b/src/physics/mhd.cpp
@@ -0,0 +1,797 @@
+#include <cmath>
+#include <vector>
+
+#include <Sacado.hpp>
+#include <deal.II/differentiation/ad/sacado_math.h>
+#include <deal.II/differentiation/ad/sacado_number_types.h>
+#include <deal.II/differentiation/ad/sacado_product_types.h>
+
+#include "physics.h"
+#include "mhd.h"
+
+
+namespace PHiLiP {
+namespace Physics {
+
+template <int dim, int nstate, typename real>
+std::array<real,nstate> MHD<dim,nstate,real>
+::source_term (
+    const dealii::Point<dim,double> &/*pos*/,
+    const std::array<real,nstate> &/*conservative_soln*/) const
+{
+    std::array<real,nstate> source_term;
+    for (int s=0; s<nstate; s++) {
+        source_term[s] = 0;
+    }
+
+    return source_term;
+}
+
+template <int dim, int nstate, typename real>
+inline std::array<real,nstate> MHD<dim,nstate,real>
+::convert_conservative_to_primitive ( const std::array<real,nstate> &conservative_soln ) const
+{
+    std::array<real, nstate> primitive_soln;
+
+    real density = conservative_soln[0];
+    dealii::Tensor<1,dim,real> vel = compute_velocities (conservative_soln);
+    real pressure = compute_pressure (conservative_soln);
+
+    primitive_soln[0] = density;
+    for (int d=0; d<dim; ++d) {
+        primitive_soln[1+d] = vel[d];
+    }
+    primitive_soln[nstate-1] = pressure;
+    return primitive_soln;
+}
+
+template <int dim, int nstate, typename real>
+inline std::array<real,nstate> MHD<dim,nstate,real>
+::convert_primitive_to_conservative ( const std::array<real,nstate> &primitive_soln ) const
+{
+
+    const real density = primitive_soln[0];
+    const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive(primitive_soln);
+
+    std::array<real, nstate> conservative_soln;
+    conservative_soln[0] = density;
+    for (int d=0; d<dim; ++d) {
+        conservative_soln[1+d] = density*velocities[d];
+    }
+    conservative_soln[nstate-1] = compute_total_energy(primitive_soln);
+
+    return conservative_soln;
+}
+
+//template <int dim, int nstate, typename real>
+//inline dealii::Tensor<1,dim,double> MHD<dim,nstate,real>::compute_velocities_inf() const
+//{
+//    dealii::Tensor<1,dim,double> velocities;
+//    return velocities;
+//}
+
+template <int dim, int nstate, typename real>
+inline dealii::Tensor<1,dim,real> MHD<dim,nstate,real>
+::compute_velocities ( const std::array<real,nstate> &conservative_soln ) const
+{
+    const real density = conservative_soln[0];
+    dealii::Tensor<1,dim,real> vel;
+    for (int d=0; d<dim; ++d) { vel[d] = conservative_soln[1+d]/density; }
+    return vel;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_velocity_squared ( const dealii::Tensor<1,dim,real> &velocities ) const
+{
+    real vel2 = 0.0;
+    for (int d=0; d<dim; d++) { vel2 = vel2 + velocities[d]*velocities[d]; }
+    return vel2;
+}
+
+template <int dim, int nstate, typename real>
+inline dealii::Tensor<1,dim,real> MHD<dim,nstate,real>
+::extract_velocities_from_primitive ( const std::array<real,nstate> &primitive_soln ) const
+{
+    dealii::Tensor<1,dim,real> velocities;
+    for (int d=0; d<dim; d++) { velocities[d] = primitive_soln[1+d]; }
+    return velocities;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_total_energy ( const std::array<real,nstate> &primitive_soln ) const
+{
+    const real density = primitive_soln[0];
+    const real pressure = primitive_soln[nstate-1];
+    const dealii::Tensor<1,dim,real> velocities = extract_velocities_from_primitive(primitive_soln);
+    const real vel2 = compute_velocity_squared(velocities);
+
+    const real tot_energy = pressure / gamm1 + 0.5*density*vel2;
+    return tot_energy;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_entropy_measure ( const std::array<real,nstate> &conservative_soln ) const
+{
+    const real density = conservative_soln[0];
+    const real pressure = compute_pressure(conservative_soln);
+    const real entropy_measure = pressure*pow(density,-gam);
+    return entropy_measure;
+}
+
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_specific_enthalpy ( const std::array<real,nstate> &conservative_soln, const real pressure ) const
+{
+    const real density = conservative_soln[0];
+    const real total_energy = conservative_soln[nstate-1];
+    const real specific_enthalpy = (total_energy+pressure)/density;
+    return specific_enthalpy;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_dimensional_temperature ( const std::array<real,nstate> &primitive_soln ) const
+{
+    const real density = primitive_soln[0];
+    const real pressure = primitive_soln[nstate-1];
+    const real temperature = gam*pressure/density;
+    return temperature;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_temperature ( const std::array<real,nstate> &primitive_soln ) const
+{
+    const real dimensional_temperature = compute_dimensional_temperature(primitive_soln);
+    const real temperature = dimensional_temperature * mach_inf_sqr;
+    return temperature;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_density_from_pressure_temperature ( const real pressure, const real temperature ) const
+{
+    const real density = gam*pressure/temperature * mach_inf_sqr;
+    return density;
+}
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_temperature_from_density_pressure ( const real density, const real pressure ) const
+{
+    const real temperature = gam*pressure/density * mach_inf_sqr;
+    return temperature;
+}
+
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_pressure ( const std::array<real,nstate> &conservative_soln ) const
+{
+    const real density = conservative_soln[0];
+    //std::cout << "density " << density << std::endl;
+
+    const real tot_energy  = conservative_soln[nstate-1];
+  //  std::cout << "tot_energy " << tot_energy << std::endl;
+
+    const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+   // std::cout << "vel1 " << vel[0] << std::endl
+   // 					 << vel[1] << std::endl
+//						 << vel[2] <<std::endl;
+
+    const real vel2 = compute_velocity_squared(vel);
+    //std::cout << "vel ^2 " << vel2 <<std::endl;
+    real pressure = gamm1*(tot_energy - 0.5*density*vel2);
+    //std::cout << "calculated pressure is" << pressure << std::endl;
+    if(pressure<0.0) {
+        std::cout<<"Cannot compute pressure..."<<std::endl;
+        std::cout<<"density "<<density<<std::endl;
+        for(int d=0;d<dim;d++) std::cout<<"vel"<<d<<" "<<vel[d]<<std::endl;
+        std::cout<<"energy "<<tot_energy<<std::endl;
+    }
+    assert(pressure>0.0);
+    //if(pressure<1e-4) pressure = 0.01;
+    return pressure;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_sound ( const std::array<real,nstate> &conservative_soln ) const
+{
+    real density = conservative_soln[0];
+    //if(density<1e-4) density = 0.01;
+    if(density<0.0) {
+        std::cout<<"density"<<density<<std::endl;
+        std::abort();
+    }
+    assert(density > 0);
+    const real pressure = compute_pressure(conservative_soln);
+    //std::cout << "pressure is" << pressure << std::endl;
+    const real sound = std::sqrt(pressure*gam/density);
+    //std::cout << "sound is " << sound << std::endl;
+    return sound;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_sound ( const real density, const real pressure ) const
+{
+    assert(density > 0);
+    const real sound = std::sqrt(pressure*gam/density);
+    return sound;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_mach_number ( const std::array<real,nstate> &conservative_soln ) const
+{
+    const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+    const real velocity = sqrt(compute_velocity_squared(vel));
+    const real sound = compute_sound (conservative_soln);
+    const real mach_number = velocity/sound;
+    return mach_number;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>
+::compute_magnetic_energy (const std::array<real,nstate> &conservative_soln) const
+{
+	const real magnetic_energy = 0;
+	for (int i = 1; i <= 3; ++i)
+		magnetic_energy += 1./2. * (conservative_soln[nstate - i] * conservative_soln[nstate - i] );
+}
+
+
+// Split form functions:
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>::
+compute_mean_density(const std::array<real,nstate> &soln_const,
+						  const std::array<real,nstate> &soln_loop) const
+{
+	return (soln_const[0] + soln_loop[0])/2.;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>::
+compute_mean_pressure(const std::array<real,nstate> &soln_const,
+    				  const std::array<real,nstate> &soln_loop) const
+{
+	real pressure_const = compute_pressure(soln_const);
+	real pressure_loop = compute_pressure(soln_loop);
+	return (pressure_const + pressure_loop)/2.;
+}
+
+template <int dim, int nstate, typename real>
+inline dealii::Tensor<1,dim,real> MHD<dim,nstate,real>::
+compute_mean_velocities(const std::array<real,nstate> &soln_const,
+    				    const std::array<real,nstate> &soln_loop) const
+{
+	dealii::Tensor<1,dim,real> vel_const = compute_velocities(soln_const);
+	dealii::Tensor<1,dim,real> vel_loop = compute_velocities(soln_loop);
+	return (vel_const + vel_loop)/2.;
+}
+
+template <int dim, int nstate, typename real>
+inline real MHD<dim,nstate,real>::
+compute_mean_specific_energy(const std::array<real,nstate> &soln_const,
+    						 const std::array<real,nstate> &soln_loop) const
+{
+	return ((soln_const[nstate-1]/soln_const[0]) + (soln_loop[nstate-1]/soln_loop[0]))/2.;
+}
+
+
+template <int dim, int nstate, typename real>
+std::array<dealii::Tensor<1,dim,real>,nstate> MHD<dim,nstate,real>
+::convective_flux (const std::array<real,nstate> &conservative_soln) const
+{
+    std::array<dealii::Tensor<1,dim,real>,nstate> conv_flux;
+    const real density = conservative_soln[0];
+    const real pressure = compute_pressure (conservative_soln);
+    const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+    const real specific_total_energy = conservative_soln[nstate-1]/conservative_soln[0];
+    const real specific_total_enthalpy = specific_total_energy + pressure/density;
+    const real magnetic_energy = compute_magnetic_energy(conservative_soln);
+
+    for (int flux_dim=0; flux_dim<dim; ++flux_dim) {
+        // Density equation
+        conv_flux[0][flux_dim] = conservative_soln[1+flux_dim];
+        // Momentum equation
+        for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
+            conv_flux[1+velocity_dim][flux_dim] = density*vel[flux_dim]*vel[velocity_dim];
+        }
+        conv_flux[1+flux_dim][flux_dim] += pressure + magnetic_energy; // Add diagonal of pressure and magnetic energy
+        // Energy equation
+        conv_flux[nstate-4][flux_dim] = density*vel[flux_dim]*specific_total_enthalpy;
+    }
+    return conv_flux;
+}
+
+template <int dim, int nstate, typename real>
+std::array<dealii::Tensor<1,dim,real>,nstate> MHD<dim, nstate, real>
+::convective_numerical_split_flux(const std::array<real,nstate> &soln_const,
+								  const std::array<real,nstate> &soln_loop) const
+{
+	std::array<dealii::Tensor<1,dim,real>,nstate> conv_num_split_flux;
+	const real mean_density = compute_mean_density(soln_const, soln_loop);
+	const real mean_pressure = compute_mean_pressure(soln_const, soln_loop);
+	const dealii::Tensor<1,dim,real> mean_velocities = compute_mean_velocities(soln_const,soln_loop);
+	const real mean_specific_energy = compute_mean_specific_energy(soln_const, soln_loop);
+
+	for (int flux_dim = 0; flux_dim < dim; ++flux_dim)
+	{
+		// Density equation
+		conv_num_split_flux[0][flux_dim] = mean_density * mean_velocities[flux_dim];//conservative_soln[1+flux_dim];
+		// Momentum equation
+		for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
+			conv_num_split_flux[1+velocity_dim][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_velocities[velocity_dim];
+		}
+		conv_num_split_flux[1+flux_dim][flux_dim] += mean_pressure; // Add diagonal of pressure
+		// Energy equation
+		conv_num_split_flux[nstate-1][flux_dim] = mean_density*mean_velocities[flux_dim]*mean_specific_energy + mean_pressure * mean_velocities[flux_dim];
+	}
+
+	return conv_num_split_flux;
+}
+
+template <int dim, int nstate, typename real>
+std::array<real,nstate> MHD<dim,nstate,real>
+::convective_normal_flux (const std::array<real,nstate> &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const
+{
+    std::array<real, nstate> conv_normal_flux;
+    const real density = conservative_soln[0];
+    const real pressure = compute_pressure (conservative_soln);
+    const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+    const real normal_vel = vel*normal;
+    const real total_energy = conservative_soln[nstate-1];
+    const real specific_total_enthalpy = (total_energy + pressure) / density;
+
+    const real rhoV = density*normal_vel;
+    // Density equation
+    conv_normal_flux[0] = rhoV;
+    // Momentum equation
+    for (int velocity_dim=0; velocity_dim<dim; ++velocity_dim){
+        conv_normal_flux[1+velocity_dim] = rhoV*vel[velocity_dim] + normal[velocity_dim] * pressure;
+    }
+    // Energy equation
+    conv_normal_flux[nstate-1] = rhoV*specific_total_enthalpy;
+    return conv_normal_flux;
+}
+
+template <int dim, int nstate, typename real>
+dealii::Tensor<2,nstate,real> MHD<dim,nstate,real>
+::convective_flux_directional_jacobian (
+    const std::array<real,nstate> &conservative_soln,
+    const dealii::Tensor<1,dim,real> &normal) const
+{
+    // See Blazek Appendix A.9 p. 429-430
+    const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+    real vel_normal = 0.0;
+    for (int d=0;d<dim;d++) { vel_normal += vel[d] * normal[d]; }
+
+    const real vel2 = compute_velocity_squared(vel);
+    const real phi = 0.5*gamm1 * vel2;
+
+    const real density = conservative_soln[0];
+    const real tot_energy = conservative_soln[nstate-1];
+    const real E = tot_energy / density;
+    const real a1 = gam*E-phi;
+    const real a2 = gam-1.0;
+    const real a3 = gam-2.0;
+
+    dealii::Tensor<2,nstate,real> jacobian;
+    for (int d=0; d<dim; ++d) {
+        jacobian[0][1+d] = normal[d];
+    }
+    for (int row_dim=0; row_dim<dim; ++row_dim) {
+        jacobian[1+row_dim][0] = normal[row_dim]*phi - vel[row_dim] * vel_normal;
+        for (int col_dim=0; col_dim<dim; ++col_dim){
+            if (row_dim == col_dim) {
+                jacobian[1+row_dim][1+col_dim] = vel_normal - a3*normal[row_dim]*vel[row_dim];
+            } else {
+                jacobian[1+row_dim][1+col_dim] = normal[col_dim]*vel[row_dim] - a2*normal[row_dim]*vel[col_dim];
+            }
+        }
+        jacobian[1+row_dim][nstate-1] = normal[row_dim]*a2;
+    }
+    jacobian[nstate-1][0] = vel_normal*(phi-a1);
+    for (int d=0; d<dim; ++d){
+        jacobian[nstate-1][1+d] = normal[d]*a1 - a2*vel[d]*vel_normal;
+    }
+    jacobian[nstate-1][nstate-1] = gam*vel_normal;
+
+    return jacobian;
+}
+
+template <int dim, int nstate, typename real>
+std::array<real,nstate> MHD<dim,nstate,real>
+::convective_eigenvalues (
+    const std::array<real,nstate> &conservative_soln,
+    const dealii::Tensor<1,dim,real> &normal) const
+{
+    const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+    std::array<real,nstate> eig;
+    real vel_dot_n = 0.0;
+    for (int d=0;d<dim;++d) { vel_dot_n += vel[d]*normal[d]; };
+    for (int i=0; i<nstate; i++) {
+        eig[i] = vel_dot_n;
+        //eig[i] = advection_speed*normal;
+
+        //eig[i] = 1.0;
+        //eig[i] = -1.0;
+    }
+    return eig;
+}
+template <int dim, int nstate, typename real>
+real MHD<dim,nstate,real>
+::max_convective_eigenvalue (const std::array<real,nstate> &conservative_soln) const
+{
+	//std::cout << "going to calculate max eig" << std::endl;
+    const dealii::Tensor<1,dim,real> vel = compute_velocities(conservative_soln);
+	//std::cout << "velocities calculated" << std::endl;
+
+    const real sound = compute_sound (conservative_soln);
+	//std::cout << "sound calculated" << std::endl;
+
+    /*const*/ real vel2 = compute_velocity_squared(vel);
+	//std::cout << "vel2 calculated" << std::endl;
+
+	if (vel2 < 0.0001)
+		vel2 = 0.0001;
+
+    const real max_eig = sqrt(vel2) + sound;
+	//std::cout << "max eig calculated" << std::endl;
+
+    return max_eig;
+}
+
+
+template <int dim, int nstate, typename real>
+std::array<dealii::Tensor<1,dim,real>,nstate> MHD<dim,nstate,real>
+::dissipative_flux (
+    const std::array<real,nstate> &/*conservative_soln*/,
+    const std::array<dealii::Tensor<1,dim,real>,nstate> &/*solution_gradient*/) const
+{
+    std::array<dealii::Tensor<1,dim,real>,nstate> diss_flux;
+    // No dissipation
+    for (int i=0; i<nstate; i++) {
+        diss_flux[i] = 0;
+    }
+    return diss_flux;
+}
+
+template <int dim, int nstate, typename real>
+void MHD<dim,nstate,real>
+::boundary_face_values (
+   const int boundary_type,
+   const dealii::Point<dim, double> &pos,
+   const dealii::Tensor<1,dim,real> &normal_int,
+   const std::array<real,nstate> &soln_int,
+   const std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_int,
+   std::array<real,nstate> &soln_bc,
+   std::array<dealii::Tensor<1,dim,real>,nstate> &soln_grad_bc) const
+{
+    // NEED TO PROVIDE AS INPUT **************************************
+    const real total_inlet_pressure = pressure_inf*pow(1.0+0.5*gamm1*mach_inf_sqr, gam/gamm1);
+    const real total_inlet_temperature = temperature_inf*pow(total_inlet_pressure/pressure_inf, gamm1/gam);
+
+    if (boundary_type == 1000) {
+        // Manufactured solution
+        std::array<real,nstate> conservative_boundary_values;
+        std::array<dealii::Tensor<1,dim,real>,nstate> boundary_gradients;
+        for (int s=0; s<nstate; s++) {
+            conservative_boundary_values[s] = this->manufactured_solution_function.value (pos, s);
+            boundary_gradients[s] = this->manufactured_solution_function.gradient (pos, s);
+        }
+        std::array<real,nstate> primitive_boundary_values = convert_conservative_to_primitive(conservative_boundary_values);
+        for (int istate=0; istate<nstate; ++istate) {
+
+            std::array<real,nstate> characteristic_dot_n = convective_eigenvalues(conservative_boundary_values, normal_int);
+            const bool inflow = (characteristic_dot_n[istate] <= 0.);
+
+            if (inflow) { // Dirichlet boundary condition
+
+                soln_bc[istate] = conservative_boundary_values[istate];
+                soln_grad_bc[istate] = soln_grad_int[istate];
+
+                // Only set the pressure and velocity
+                // primitive_boundary_values[0] = soln_int[0];;
+                // for(int d=0;d<dim;d++){
+                //    primitive_boundary_values[1+d] = soln_int[1+d]/soln_int[0];;
+                //}
+                conservative_boundary_values = convert_primitive_to_conservative(primitive_boundary_values);
+                //conservative_boundary_values[nstate-1] = soln_int[nstate-1];
+                soln_bc[istate] = conservative_boundary_values[istate];
+
+            } else { // Neumann boundary condition
+                // soln_bc[istate] = -soln_int[istate]+2*conservative_boundary_values[istate];
+                soln_bc[istate] = soln_int[istate];
+
+                // **************************************************************************************************************
+                // Note I don't know how to properly impose the soln_grad_bc to obtain an adjoint consistent scheme
+                // Currently, Neumann boundary conditions are only imposed for the linear advection
+                // Therefore, soln_grad_bc does not affect the solution
+                // **************************************************************************************************************
+                soln_grad_bc[istate] = soln_grad_int[istate];
+                //soln_grad_bc[istate] = boundary_gradients[istate];
+                //soln_grad_bc[istate] = -soln_grad_int[istate]+2*boundary_gradients[istate];
+            }
+
+            // HARDCODE DIRICHLET BC
+            soln_bc[istate] = conservative_boundary_values[istate];
+
+        }
+    } else if (boundary_type == 1001) {
+        // No penetration,
+        // Given by Algorithm II of the following paper
+        // Krivodonova, L., and Berger, M.,
+        // “High-order accurate implementation of solid wall boundary conditions in curved geometries,”
+        // Journal of Computational Physics, vol. 211, 2006, pp. 492–512.
+        const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive(soln_int);
+
+        // Copy density and pressure
+        std::array<real,nstate> primitive_boundary_values;
+        primitive_boundary_values[0] = primitive_interior_values[0];
+        primitive_boundary_values[nstate-1] = primitive_interior_values[nstate-1];
+
+        const dealii::Tensor<1,dim,real> surface_normal = -normal_int;
+        const dealii::Tensor<1,dim,real> velocities_int = extract_velocities_from_primitive(primitive_interior_values);
+        const dealii::Tensor<1,dim,real> velocities_bc = velocities_int - 2.0*(velocities_int*surface_normal)*surface_normal;
+        for (int d=0; d<dim; ++d) {
+            primitive_boundary_values[1+d] = velocities_bc[d];
+        }
+
+        soln_bc = convert_primitive_to_conservative(primitive_boundary_values);
+
+    } else if (boundary_type == 1002) {
+        // Pressure Outflow Boundary Condition (back pressure)
+        // Carlson 2011, sec. 2.4
+
+        const real back_pressure = 0.99; // Make it as an input later on
+
+        const real mach_int = compute_mach_number(soln_int);
+        const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive(soln_int);
+        const real pressure_int = primitive_interior_values[nstate-1];
+
+        const real radicant = 1.0+0.5*gamm1*mach_inf_sqr;
+        const real pressure_inlet = total_inlet_pressure * pow(radicant, -gam/gamm1);
+        const real pressure_bc = (mach_int >= 1) ? pressure_int : back_pressure*pressure_inlet;
+        const real temperature_int = compute_temperature(primitive_interior_values);
+
+        // Assign primitive boundary values
+        std::array<real,nstate> primitive_boundary_values;
+        primitive_boundary_values[0] = compute_density_from_pressure_temperature(pressure_bc, temperature_int);
+        for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = primitive_interior_values[1+d]; }
+        primitive_boundary_values[nstate-1] = pressure_bc;
+
+        soln_bc = convert_primitive_to_conservative(primitive_boundary_values);
+
+        // Supersonic, simply extrapolate
+        if (mach_int > 1.0) {
+            soln_bc = soln_int;
+        }
+
+    } else if (boundary_type == 1003) {
+        // Inflow
+        // Carlson 2011, sec. 2.2 & sec 2.9
+
+        const std::array<real,nstate> primitive_interior_values = convert_conservative_to_primitive(soln_int);
+
+        const dealii::Tensor<1,dim,real> normal = -normal_int;
+
+        const real                       density_i    = primitive_interior_values[0];
+        const dealii::Tensor<1,dim,real> velocities_i = extract_velocities_from_primitive(primitive_interior_values);
+        const real                       pressure_i   = primitive_interior_values[nstate-1];
+
+        const real                       normal_vel_i = velocities_i*normal;
+        const real                       sound_i      = compute_sound(soln_int);
+        //const real                       mach_i       = std::abs(normal_vel_i)/sound_i;
+
+        //const dealii::Tensor<1,dim,real> velocities_o = velocities_inf;
+        //const real                       normal_vel_o = velocities_o*normal;
+        //const real                       sound_o      = sound_inf;
+        //const real                       mach_o       = mach_inf;
+
+        if(mach_inf < 1.0) {
+            //std::cout << "Subsonic inflow, mach=" << mach_i << std::endl;
+            // Subsonic inflow, sec 2.7
+
+            // Want to solve for c_b (sound_bc), to then solve for U (velocity_magnitude_bc) and M_b (mach_bc)
+            // Eq. 37
+            const real riemann_pos = normal_vel_i + 2.0*sound_i/gamm1;
+            // Could evaluate enthalpy from primitive like eq.36, but easier to use the following
+            const real specific_total_energy = soln_int[nstate-1]/density_i;
+            const real specific_total_enthalpy = specific_total_energy + pressure_i/density_i;
+            // Eq. 43
+            const real a = 1.0+2.0/gamm1;
+            const real b = -2.0*riemann_pos;
+            const real c = 0.5*gamm1 * (riemann_pos*riemann_pos - 2.0*specific_total_enthalpy);
+            // Eq. 42
+            const real term1 = -0.5*b/a;
+            const real term2= 0.5*sqrt(b*b-4.0*a*c)/a;
+            const real sound_bc1 = term1 + term2;
+            const real sound_bc2 = term1 - term2;
+            // Eq. 44
+            const real sound_bc  = std::max(sound_bc1, sound_bc2);
+            // Eq. 45
+            //const real velocity_magnitude_bc = 2.0*sound_bc/gamm1 - riemann_pos;
+            const real velocity_magnitude_bc = riemann_pos - 2.0*sound_bc/gamm1;
+            const real mach_bc = velocity_magnitude_bc/sound_bc;
+            // Eq. 46
+            const real radicant = 1.0+0.5*gamm1*mach_bc*mach_bc;
+            const real pressure_bc = total_inlet_pressure * pow(radicant, -gam/gamm1);
+            const real temperature_bc = total_inlet_temperature * pow(radicant, -1.0);
+            //std::cout << " pressure_bc " << pressure_bc << "pressure_inf" << pressure_inf << std::endl;
+            //std::cout << " temperature_bc " << temperature_bc << "temperature_inf" << temperature_inf << std::endl;
+            //
+
+            const real density_bc  = compute_density_from_pressure_temperature(pressure_bc, temperature_bc);
+            std::array<real,nstate> primitive_boundary_values;
+            primitive_boundary_values[0] = density_bc;
+            for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = velocity_magnitude_bc*normal[d]; }
+            primitive_boundary_values[nstate-1] = pressure_bc;
+            soln_bc = convert_primitive_to_conservative(primitive_boundary_values);
+
+            //std::cout << " entropy_bc " << compute_entropy_measure(soln_bc) << "entropy_inf" << entropy_inf << std::endl;
+
+        } else {
+            // Supersonic inflow, sec 2.9
+            // Specify all quantities through
+            // total_inlet_pressure, total_inlet_temperature, mach_inf & angle_of_attack
+            //std::cout << "Supersonic inflow, mach=" << mach_i << std::endl;
+            const real radicant = 1.0+0.5*gamm1*mach_inf_sqr;
+            const real static_inlet_pressure    = total_inlet_pressure * pow(radicant, -gam/gamm1);
+            const real static_inlet_temperature = total_inlet_temperature * pow(radicant, -1.0);
+
+            const real pressure_bc = static_inlet_pressure;
+            const real temperature_bc = static_inlet_temperature;
+            const real density_bc  = compute_density_from_pressure_temperature(pressure_bc, temperature_bc);
+            const real sound_bc = sqrt(gam * pressure_bc / density_bc);
+            const real velocity_magnitude_bc = mach_inf * sound_bc;
+
+            // Assign primitive boundary values
+            std::array<real,nstate> primitive_boundary_values;
+            primitive_boundary_values[0] = density_bc;
+            for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = -velocity_magnitude_bc*normal_int[d]; } // minus since it's inflow
+            primitive_boundary_values[nstate-1] = pressure_bc;
+            soln_bc = convert_primitive_to_conservative(primitive_boundary_values);
+            //std::cout << "Inlet density : " << density_bc << std::endl;
+            //std::cout << "Inlet vel_x   : " << primitive_boundary_values[1] << std::endl;
+            //std::cout << "Inlet vel_y   : " << primitive_boundary_values[2] << std::endl;
+            //std::cout << "Inlet pressure: " << pressure_bc << std::endl;
+        }
+
+    } else if (boundary_type == 1004) {
+        // Farfield boundary condition
+        const real density_bc = density_inf;
+        const real pressure_bc = 1.0/(gam*mach_inf_sqr);
+        std::array<real,nstate> primitive_boundary_values;
+        primitive_boundary_values[0] = density_bc;
+        for (int d=0;d<dim;d++) { primitive_boundary_values[1+d] = velocities_inf[d]; } // minus since it's inflow
+        primitive_boundary_values[nstate-1] = pressure_bc;
+        soln_bc = convert_primitive_to_conservative(primitive_boundary_values);
+        //std::cout << "Density inf " << soln_bc[0] << std::endl;
+        //std::cout << "momxinf " << soln_bc[1] << std::endl;
+        //std::cout << "momxinf " << soln_bc[2] << std::endl;
+        //std::cout << "energy inf " << soln_bc[3] << std::endl;
+    } else{
+        std::cout << "Invalid boundary_type: " << boundary_type << std::endl;
+        std::abort();
+    }
+}
+
+template <int dim, int nstate, typename real>
+dealii::Vector<double> MHD<dim,nstate,real>::post_compute_derived_quantities_vector (
+    const dealii::Vector<double>              &uh,
+    const std::vector<dealii::Tensor<1,dim> > &duh,
+    const std::vector<dealii::Tensor<2,dim> > &dduh,
+    const dealii::Tensor<1,dim>               &normals,
+    const dealii::Point<dim>                  &evaluation_points) const
+{
+    std::vector<std::string> names = post_get_names ();
+    dealii::Vector<double> computed_quantities = PhysicsBase<dim,nstate,real>::post_compute_derived_quantities_vector ( uh, duh, dduh, normals, evaluation_points);
+    unsigned int current_data_index = computed_quantities.size() - 1;
+    computed_quantities.grow_or_shrink(names.size());
+    if constexpr (std::is_same<real,double>::value) {
+
+        std::array<double, nstate> conservative_soln;
+        for (unsigned int s=0; s<nstate; ++s) {
+            conservative_soln[s] = uh(s);
+        }
+        const std::array<double, nstate> primitive_soln = convert_conservative_to_primitive(conservative_soln);
+
+        // Density
+        computed_quantities(++current_data_index) = primitive_soln[0];
+        // Velocities
+        for (unsigned int d=0; d<dim; ++d) {
+            computed_quantities(++current_data_index) = primitive_soln[1+d];
+        }
+        // Momentum
+        for (unsigned int d=0; d<dim; ++d) {
+            computed_quantities(++current_data_index) = conservative_soln[1+d];
+        }
+        // Energy
+        computed_quantities(++current_data_index) = conservative_soln[nstate-1];
+        // Pressure
+        computed_quantities(++current_data_index) = primitive_soln[nstate-1];
+        // Pressure
+        computed_quantities(++current_data_index) = compute_temperature(primitive_soln);
+        // Entropy generation
+        computed_quantities(++current_data_index) = compute_entropy_measure(conservative_soln) - entropy_inf;
+        // Mach Number
+        computed_quantities(++current_data_index) = compute_mach_number(conservative_soln);
+
+    }
+    if (computed_quantities.size()-1 != current_data_index) {
+        std::cout << " Did not assign a value to all the data. Missing " << computed_quantities.size() - current_data_index << " variables."
+                  << " If you added a new output variable, make sure the names and DataComponentInterpretation match the above. "
+                  << std::endl;
+    }
+
+    return computed_quantities;
+}
+
+template <int dim, int nstate, typename real>
+std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> MHD<dim,nstate,real>
+::post_get_data_component_interpretation () const
+{
+    namespace DCI = dealii::DataComponentInterpretation;
+    std::vector<DCI::DataComponentInterpretation> interpretation = PhysicsBase<dim,nstate,real>::post_get_data_component_interpretation (); // state variables
+    interpretation.push_back (DCI::component_is_scalar); // Density
+    for (unsigned int d=0; d<dim; ++d) {
+        interpretation.push_back (DCI::component_is_part_of_vector); // Velocity
+    }
+    for (unsigned int d=0; d<dim; ++d) {
+        interpretation.push_back (DCI::component_is_part_of_vector); // Momentum
+    }
+    interpretation.push_back (DCI::component_is_scalar); // Energy
+    interpretation.push_back (DCI::component_is_scalar); // Pressure
+    interpretation.push_back (DCI::component_is_scalar); // Temperature
+    interpretation.push_back (DCI::component_is_scalar); // Entropy generation
+    interpretation.push_back (DCI::component_is_scalar); // Mach number
+
+    std::vector<std::string> names = post_get_names();
+    if (names.size() != interpretation.size()) {
+        std::cout << "Number of DataComponentInterpretation is not the same as number of names for output file" << std::endl;
+    }
+    return interpretation;
+}
+
+
+template <int dim, int nstate, typename real>
+std::vector<std::string> MHD<dim,nstate,real> ::post_get_names () const
+{
+    std::vector<std::string> names = PhysicsBase<dim,nstate,real>::post_get_names ();
+    names.push_back ("density");
+    for (unsigned int d=0; d<dim; ++d) {
+      names.push_back ("velocity");
+    }
+    for (unsigned int d=0; d<dim; ++d) {
+      names.push_back ("momentum");
+    }
+    names.push_back ("energy");
+    names.push_back ("pressure");
+    names.push_back ("temperature");
+
+    names.push_back ("entropy_generation");
+    names.push_back ("mach_number");
+    return names;
+}
+
+template <int dim, int nstate, typename real>
+dealii::UpdateFlags MHD<dim,nstate,real>
+::post_get_needed_update_flags () const
+{
+    //return update_values | update_gradients;
+    return dealii::update_values;
+}
+
+// Instantiate explicitly
+template class MHD < PHILIP_DIM, 2*PHILIP_DIM+2, double >;
+template class MHD < PHILIP_DIM, 2*PHILIP_DIM+2, Sacado::Fad::DFad<double>  >;
+
+} // Physics namespace
+} // PHiLiP namespace
diff --git a/src/physics/mhd.h b/src/physics/mhd.h
new file mode 100644
index 000000000..2df54a10d
--- /dev/null
+++ b/src/physics/mhd.h
@@ -0,0 +1,237 @@
+#ifndef __MHD__
+#define __MHD__
+
+#include <deal.II/base/tensor.h>
+#include "physics.h"
+
+namespace PHiLiP {
+namespace Physics {
+
+/// Magnetohydrodynamics (MHD) equations. Derived from PhysicsBase
+/** Only 2D and 3D
+ *  State variable and convective fluxes given by
+ *
+ *  \f[
+ *  \mathbf{w} =
+ *  \begin{bmatrix} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \\ \rho E \end{bmatrix}
+ *  , \qquad
+ *  \mathbf{F}_{conv} =
+ *  \begin{bmatrix}
+ *      \mathbf{f}^x_{conv}, \mathbf{f}^y_{conv}, \mathbf{f}^z_{conv}
+ *  \end{bmatrix}
+ *  =
+ *  \begin{bmatrix}
+ *  \begin{bmatrix}
+ *  \rho v_1 \\
+ *  \rho v_1 v_1 + p \\
+ *  \rho v_1 v_2     \\
+ *  \rho v_1 v_3     \\
+ *  v_1 (\rho e+p)
+ *  \end{bmatrix}
+ *  ,
+ *  \begin{bmatrix}
+ *  \rho v_2 \\
+ *  \rho v_1 v_2     \\
+ *  \rho v_2 v_2 + p \\
+ *  \rho v_2 v_3     \\
+ *  v_2 (\rho e+p)
+ *  \end{bmatrix}
+ *  ,
+ *  \begin{bmatrix}
+ *  \rho v_3 \\
+ *  \rho v_1 v_3     \\
+ *  \rho v_2 v_3     \\
+ *  \rho v_3 v_3 + p \\
+ *  v_3 (\rho e+p)
+ *  \end{bmatrix}
+ *  \end{bmatrix} \f]
+ *
+ *  where, \f$ E \f$ is the specific total energy and \f$ e \f$ is the specific internal
+ *  energy, related by
+ *  \f[
+ *      E = e + |V|^2 / 2
+ *  \f]
+ *  For a calorically perfect gas
+ *
+ *  \f[
+ *  p=(\gamma -1)(\rho e-\frac{1}{2}\rho \|\mathbf{v}\|)
+ *  \f]
+ *
+ *  Dissipative flux \f$ \mathbf{F}_{diss} = \mathbf{0} \f$
+ *
+ *  Source term \f$ s(\mathbf{x}) \f$
+ *
+ *  Equation:
+ *  \f[ \boldsymbol{\nabla} \cdot
+ *         (  \mathbf{F}_{conv}( w )
+ *          + \mathbf{F}_{diss}( w, \boldsymbol{\nabla}(w) )
+ *      = s(\mathbf{x})
+ *  \f]
+ *
+ *
+ *  Still need to provide functions to un-non-dimensionalize the variables.
+ *  Like, given density_inf
+ */
+template <int dim, int nstate, typename real>
+class MHD : public PhysicsBase <dim, nstate, real>
+{
+public:
+    /// Constructor
+    MHD (const double gamma_gas)
+    : gam(gamma_gas)
+    , gamm1(gam-1.0)
+    {
+        static_assert(nstate==2*dim+2, "Physics::MHD() should be created with nstate=dim+2");
+
+    };
+    /// Destructor
+    ~MHD ()
+    {};
+
+    /// Constant heat capacity ratio of air
+    const double gam;
+    /// Gamma-1.0 used often
+    const double gamm1;
+
+    std::array<real,nstate> manufactured_solution (const dealii::Point<dim,double> &pos) const;
+
+    /// Convective flux: \f$ \mathbf{F}_{conv} \f$
+    std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (
+        const std::array<real,nstate> &conservative_soln) const;
+
+    std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux (
+           		const std::array<real,nstate> &soln_const, const std::array<real,nstate> &soln_loop) const;
+
+
+    std::array<real,nstate> convective_normal_flux (const std::array<real,nstate> &conservative_soln, const dealii::Tensor<1,dim,real> &normal) const;
+
+    /// Convective flux Jacobian: \f$ \frac{\partial \mathbf{F}_{conv}}{\partial w} \cdot \mathbf{n} \f$
+    dealii::Tensor<2,nstate,real> convective_flux_directional_jacobian (
+        const std::array<real,nstate> &conservative_soln,
+        const dealii::Tensor<1,dim,real> &normal) const;
+
+    /// Spectral radius of convective term Jacobian is 'c'
+    std::array<real,nstate> convective_eigenvalues (
+        const std::array<real,nstate> &/*conservative_soln*/,
+        const dealii::Tensor<1,dim,real> &/*normal*/) const;
+
+    /// Maximum convective eigenvalue used in Lax-Friedrichs
+    real max_convective_eigenvalue (const std::array<real,nstate> &soln) const;
+
+    /// Dissipative flux: 0
+    std::array<dealii::Tensor<1,dim,real>,nstate> dissipative_flux (
+        const std::array<real,nstate> &conservative_soln,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &solution_gradient) const;
+
+    /// Source term is zero or depends on manufactured solution
+    std::array<real,nstate> source_term (
+        const dealii::Point<dim,double> &pos,
+        const std::array<real,nstate> &conservative_soln) const;
+
+    /// Given conservative variables [density, [momentum], total energy],
+    /// returns primitive variables [density, [velocities], pressure].
+    ///
+    /// Opposite of convert_primitive_to_conservative
+    std::array<real,nstate> convert_conservative_to_primitive ( const std::array<real,nstate> &conservative_soln ) const;
+
+    /// Given primitive variables [density, [velocities], pressure],
+    /// returns conservative variables [density, [momentum], total energy].
+    ///
+    /// Opposite of convert_primitive_to_conservative
+    std::array<real,nstate> convert_primitive_to_conservative ( const std::array<real,nstate> &primitive_soln ) const;
+
+    /// Evaluate pressure from conservative variables
+    real compute_pressure ( const std::array<real,nstate> &conservative_soln ) const;
+
+    /// Evaluate Magnetic Energy
+    real compute_magnetic_energy (const std::array<real,nstate> &conservative_soln) const;
+
+    /// Evaluate pressure from conservative variables
+    real compute_pressure_from_enthalpy ( const std::array<real,nstate> &conservative_soln ) const;
+
+    /// Evaluate pressure from conservative variables
+    real compute_specific_enthalpy ( const std::array<real,nstate> &conservative_soln, const real pressure) const;
+
+    /// Evaluate speed of sound from conservative variables
+    real compute_sound ( const std::array<real,nstate> &conservative_soln ) const;
+    /// Evaluate speed of sound from density and pressure
+    real compute_sound ( const real density, const real pressure ) const;
+
+    /// Evaluate velocities from conservative variables
+    dealii::Tensor<1,dim,real> compute_velocities ( const std::array<real,nstate> &conservative_soln ) const;
+    /// Given the velocity vector \f$ \mathbf{u} \f$, returns the dot-product  \f$ \mathbf{u} \cdot \mathbf{u} \f$
+    real compute_velocity_squared ( const dealii::Tensor<1,dim,real> &velocities ) const;
+
+    /// Given primitive variables, returns velocities.
+    dealii::Tensor<1,dim,real> extract_velocities_from_primitive ( const std::array<real,nstate> &primitive_soln ) const;
+    /// Given primitive variables, returns total energy
+    real compute_total_energy ( const std::array<real,nstate> &primitive_soln ) const;
+
+    /// Evaluate entropy from conservative variables
+    /** Note that it is not the actual entropy since it's missing some constants.
+     *  Used to check entropy convergence
+     *  See discussion in
+     *  https://physics.stackexchange.com/questions/116779/entropy-is-constant-how-to-express-this-equation-in-terms-of-pressure-and-densi?answertab=votes#tab-top
+     */
+    real compute_entropy_measure ( const std::array<real,nstate> &conservative_soln ) const;
+
+    /// Given conservative variables, returns Mach number
+    real compute_mach_number ( const std::array<real,nstate> &conservative_soln ) const;
+
+    /// Given primitive variables, returns DIMENSIONALIZED temperature using the equation of state
+    real compute_dimensional_temperature ( const std::array<real,nstate> &primitive_soln ) const;
+
+    /// Given primitive variables, returns NON-DIMENSIONALIZED temperature using free-stream non-dimensionalization
+    /** See the book I do like CFD, sec 4.14.2 */
+    real compute_temperature ( const std::array<real,nstate> &primitive_soln ) const;
+
+    /// Given pressure and temperature, returns NON-DIMENSIONALIZED density using free-stream non-dimensionalization
+    /** See the book I do like CFD, sec 4.14.2 */
+    real compute_density_from_pressure_temperature ( const real pressure, const real temperature ) const;
+
+    /// Given density and pressure, returns NON-DIMENSIONALIZED temperature using free-stream non-dimensionalization
+    /** See the book I do like CFD, sec 4.14.2 */
+    real compute_temperature_from_density_pressure ( const real density, const real pressure ) const;
+
+    ///These functions are only relevant to the split form. The Euler split form is that of Kennedy & Gruber.
+    /** Refer to Gassner's paper for more information:  */
+    real compute_mean_density(const std::array<real,nstate> &soln_const,
+    						  const std::array<real,nstate> &soln_loop) const;
+
+    real compute_mean_pressure(const std::array<real,nstate> &soln_const,
+        					   const std::array<real,nstate> &soln_loop) const;
+
+    dealii::Tensor<1,dim,real> compute_mean_velocities(const std::array<real,nstate> &soln_const,
+        						  	  	  	  	  	   const std::array<real,nstate> &soln_loop) const;
+
+    real compute_mean_specific_energy(const std::array<real,nstate> &soln_const,
+        						  const std::array<real,nstate> &soln_loop) const;
+
+    void boundary_face_values (
+        const int /*boundary_type*/,
+        const dealii::Point<dim, double> &/*pos*/,
+        const dealii::Tensor<1,dim,real> &/*normal*/,
+        const std::array<real,nstate> &/*soln_int*/,
+        const std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_int*/,
+        std::array<real,nstate> &/*soln_bc*/,
+        std::array<dealii::Tensor<1,dim,real>,nstate> &/*soln_grad_bc*/) const;
+
+    virtual dealii::Vector<double> post_compute_derived_quantities_vector (
+        const dealii::Vector<double>      &uh,
+        const std::vector<dealii::Tensor<1,dim> > &duh,
+        const std::vector<dealii::Tensor<2,dim> > &dduh,
+        const dealii::Tensor<1,dim>                  &normals,
+        const dealii::Point<dim>                  &evaluation_points) const;
+    virtual std::vector<std::string> post_get_names () const;
+    virtual std::vector<dealii::DataComponentInterpretation::DataComponentInterpretation> post_get_data_component_interpretation () const;
+    virtual dealii::UpdateFlags post_get_needed_update_flags () const;
+protected:
+
+
+};
+
+} // Physics namespace
+} // PHiLiP namespace
+
+#endif
+
diff --git a/src/physics/physics.h b/src/physics/physics.h
index 90a51b470..4112d9d3d 100644
--- a/src/physics/physics.h
+++ b/src/physics/physics.h
@@ -44,6 +44,10 @@ class PhysicsBase
     virtual std::array<dealii::Tensor<1,dim,real>,nstate> convective_flux (
         const std::array<real,nstate> &solution) const = 0;
 
+    /// Convective Numerical Split Flux for split form
+    virtual std::array<dealii::Tensor<1,dim,real>,nstate> convective_numerical_split_flux (
+    		const std::array<real,nstate> &soln_const, const std::array<real,nstate> &soln_loop) const = 0;
+
     /// Spectral radius of convective term Jacobian.
     /// Used for scalar dissipation
     virtual std::array<real,nstate> convective_eigenvalues (
@@ -114,9 +118,6 @@ class PhysicsBase
      */
     dealii::Tensor<2,dim,double> diffusion_tensor;
 };
-
-
-
 } // Physics namespace
 } // PHiLiP namespace
 
diff --git a/src/physics/physics_factory.cpp b/src/physics/physics_factory.cpp
index 17dc2e6e0..665c4753b 100644
--- a/src/physics/physics_factory.cpp
+++ b/src/physics/physics_factory.cpp
@@ -5,6 +5,7 @@
 #include "convection_diffusion.h"
 #include "burgers.h"
 #include "euler.h"
+//#include "mhd.h"
 
 namespace PHiLiP {
 namespace Physics {
@@ -35,7 +36,9 @@ ::create_Physics(const Parameters::AllParameters *const parameters_input)
                                                                ,parameters_input->euler_param.side_slip_angle);
         }
 
-    }
+    } //else if (pde_type == PDE_enum::mhd) {
+    	//if constexpr (nstate == 2*dim + 2) return std::make_shared < MHD<dim,nstate,real> > ();
+   // }
     std::cout << "Can't create PhysicsBase, invalid PDE type: " << pde_type << std::endl;
     assert(0==1 && "Can't create PhysicsBase, invalid PDE type");
     return nullptr;
diff --git a/src/temporary_test.cpp b/src/temporary_test.cpp
index a48bea0fa..937463030 100644
--- a/src/temporary_test.cpp
+++ b/src/temporary_test.cpp
@@ -1,3 +1,4 @@
+// hello
 #include <Sacado.hpp>
 #include <deal.II/differentiation/ad/sacado_product_types.h>
 #include <deal.II/base/tensor.h>
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
index b7e939819..17ae85eb4 100644
--- a/src/tests/CMakeLists.txt
+++ b/src/tests/CMakeLists.txt
@@ -26,7 +26,7 @@ foreach(dim RANGE 1 3)
     # Setup target with deal.II
     DEAL_II_SETUP_TARGET(${TestsLib})
 
-    unset(Testsib)
+    unset(TestsLib)
     unset(DiscontinuousGalerkinLib)
     unset(ODESolverLib)
     unset(NumericalFluxLib)
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 88577f66e..9fca3c2b5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -9,5 +9,8 @@ add_subdirectory(diffusion_implicit)
 add_subdirectory(convection_diffusion_implicit)
 add_subdirectory(advection_vector_implicit)
 add_subdirectory(burgers_inviscid_implicit)
+add_subdirectory(burgers_split_stability)
+add_subdirectory(euler_split_inviscid_taylor_green_vortex)
+add_subdirectory(advection_explicit_periodic)
 
 add_subdirectory(euler_integration)
diff --git a/tests/advection_explicit_periodic/2D_advection_explicit_periodic.prm b/tests/advection_explicit_periodic/2D_advection_explicit_periodic.prm
new file mode 100644
index 000000000..8b4798e6d
--- /dev/null
+++ b/tests/advection_explicit_periodic/2D_advection_explicit_periodic.prm
@@ -0,0 +1,30 @@
+# -------------------
+# Number of dimensions
+set dimension = 2
+
+set use_weak_form = true
+
+set use_collocated_nodes = false
+
+set use_split_form = false
+
+set use_periodic_bc = true
+
+# The PDE we want to solve
+set pde_type = advection
+
+set conv_num_flux = lax_friedrichs
+
+subsection ODE solver
+
+  set ode_output = verbose
+  
+  set nonlinear_max_iterations = 500
+
+  set print_iteration_modulo = 100
+
+  set ode_solver_type = explicit
+
+  set initial_time_step = 0.001
+
+end
diff --git a/tests/advection_explicit_periodic/CMakeLists.txt b/tests/advection_explicit_periodic/CMakeLists.txt
new file mode 100644
index 000000000..3a46ce27e
--- /dev/null
+++ b/tests/advection_explicit_periodic/CMakeLists.txt
@@ -0,0 +1,52 @@
+set(TEST_SRC
+    advection_explicit_periodic.cpp
+    )
+
+foreach(dim RANGE 2 2)
+
+    
+    # Output executable
+    string(CONCAT TEST_TARGET ${dim}D_advection_explicit_periodic)
+    message("Adding executable " ${TEST_TARGET} " with files " ${TEST_SRC} "\n")
+    add_executable(${TEST_TARGET} ${TEST_SRC})
+    # Replace occurences of PHILIP_DIM with 1, 2, or 3 in the code
+    target_compile_definitions(${TEST_TARGET} PRIVATE PHILIP_DIM=${dim})
+
+    # Compile this executable when 'make unit_tests'
+    add_dependencies(unit_tests ${TEST_TARGET})
+    add_dependencies(${dim}D ${TEST_TARGET})
+
+    # Library dependency
+    set(ParameterLib ParametersLibrary)
+    string(CONCAT PhysicsLib Physics_${dim}D)
+    string(CONCAT NumericalFluxLib NumericalFlux_${dim}D)
+    string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) 
+    string(CONCAT ODESolverLib ODESolver_${dim}D)
+
+    target_link_libraries(${TEST_TARGET} ${ParameterLib})
+    target_link_libraries(${TEST_TARGET} ${PhysicsLib})
+    target_link_libraries(${TEST_TARGET} ${NumericalFluxLib})
+    target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib})
+    target_link_libraries(${TEST_TARGET} ${ODESolverLib})
+
+    #target_link_libraries(${TEST_TARGET} m mpi mpi_cxx)
+
+    # Setup target with deal.II
+    DEAL_II_SETUP_TARGET(${TEST_TARGET})
+    
+    configure_file(2D_advection_explicit_periodic.prm 2D_advection_explicit_periodic.prm COPYONLY)
+    add_test(
+      NAME ${TEST_TARGET}
+      COMMAND ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} -i ${CMAKE_CURRENT_BINARY_DIR}/2D_advection_explicit_periodic.prm
+      WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+    )
+
+    unset(dim)
+    unset(TEST_TARGET)
+    unset(PhysicsLib)
+    unset(NumericalFluxLib)
+    unset(ParameterLib)
+    unset(DiscontinuousGalerkinLib)
+    unset(ODESolverLib)
+
+endforeach()
diff --git a/tests/advection_explicit_periodic/advection_explicit_periodic.cpp b/tests/advection_explicit_periodic/advection_explicit_periodic.cpp
new file mode 100644
index 000000000..dc5d2657e
--- /dev/null
+++ b/tests/advection_explicit_periodic/advection_explicit_periodic.cpp
@@ -0,0 +1,203 @@
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/function.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/base/numbers.h>
+#include <deal.II/base/function_parser.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+
+#include "parameters/all_parameters.h"
+#include "parameters/parameters.h"
+#include "numerical_flux/numerical_flux.h"
+#include "physics/physics_factory.h"
+#include "physics/physics.h"
+#include "dg/dg.h"
+#include "ode_solver/ode_solver.h"
+
+#include<fenv.h>
+
+//using PDEType  = PHiLiP::Parameters::AllParameters::PartialDifferentialEquation;
+//using ConvType = PHiLiP::Parameters::AllParameters::ConvectiveNumericalFlux;
+//using DissType = PHiLiP::Parameters::AllParameters::DissipativeNumericalFlux;
+//
+//
+//const double TOLERANCE = 1E-12;
+
+
+template <int dim, int nstate>
+class AdvectionPeriodic
+{
+public:
+	AdvectionPeriodic() = delete;
+	AdvectionPeriodic(const PHiLiP::Parameters::AllParameters *const parameters_input);
+	int run_test();
+
+
+
+
+private:
+    const PHiLiP::Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters
+    const MPI_Comm mpi_communicator;
+    dealii::ConditionalOStream pcout;
+};
+
+template <int dim, int nstate>
+AdvectionPeriodic<dim, nstate>::AdvectionPeriodic(const PHiLiP::Parameters::AllParameters *const parameters_input)
+:
+all_parameters(parameters_input)
+, mpi_communicator(MPI_COMM_WORLD)
+, pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0)
+{}
+
+
+template <int dim, int nstate>
+int AdvectionPeriodic<dim, nstate>::run_test()
+{
+#if PHILIP_DIM==1 // dealii::parallel::distributed::Triangulation<dim> does not work for 1D
+            dealii::Triangulation<dim> grid(
+                typename dealii::Triangulation<dim>::MeshSmoothing(
+                    dealii::Triangulation<dim>::smoothing_on_refinement |
+                    dealii::Triangulation<dim>::smoothing_on_coarsening));
+#else
+            dealii::parallel::distributed::Triangulation<dim> grid(
+                this->mpi_communicator);
+#endif
+	//dealii::Triangulation<dim> grid;
+
+	double left = 0.0;
+	double right = 2.0;
+	const bool colorize = true;
+	int n_refinements = 5;
+	unsigned int poly_degree = 2;
+	dealii::GridGenerator::hyper_cube(grid, left, right, colorize);
+
+	std::vector<dealii::GridTools::PeriodicFacePair<typename dealii::parallel::distributed::Triangulation<PHILIP_DIM>::cell_iterator> > matched_pairs;
+		dealii::GridTools::collect_periodic_faces(grid,0,1,0,matched_pairs);
+		dealii::GridTools::collect_periodic_faces(grid,2,3,1,matched_pairs);
+		//dealii::GridTools::collect_periodic_faces(grid,4,5,2,matched_pairs);
+		grid.add_periodicity(matched_pairs);
+
+	grid.refine_global(n_refinements);
+
+	std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree);
+	dg->set_triangulation(&grid);
+	dg->allocate_system ();
+
+//	for (auto current_cell = dg->dof_handler.begin_active(); current_cell != dg->dof_handler.end(); ++current_cell) {
+//		 if (!current_cell->is_locally_owned()) continue;
+//
+//		 dg->fe_values_volume.reinit(current_cell);
+//		 int cell_index = current_cell->index();
+//		 std::cout << "cell number " << cell_index << std::endl;
+//		 for (unsigned int face_no = 0; face_no < dealii::GeometryInfo<PHILIP_DIM>::faces_per_cell; ++face_no)
+//		 {
+//			 if (current_cell->face(face_no)->at_boundary())
+//		     {
+//				 std::cout << "face " << face_no << " is at boundary" << std::endl;
+//		         typename dealii::DoFHandler<PHILIP_DIM>::active_cell_iterator neighbor = current_cell->neighbor_or_periodic_neighbor(face_no);
+//		         std::cout << "the neighbor is " << neighbor->index() << std::endl;
+//		     }
+//		 }
+//
+//	}
+
+	std::cout << "Implement initial conditions" << std::endl;
+	dealii::FunctionParser<2> initial_condition;
+	std::string variables = "x,y";
+	std::map<std::string,double> constants;
+	constants["pi"] = dealii::numbers::PI;
+	std::string expression = "exp( -( 20*(x-1)*(x-1) + 20*(y-1)*(y-1) ) )";//"sin(pi*x)*sin(pi*y)";
+	initial_condition.initialize(variables,
+	                             expression,
+	                             constants);
+	dealii::VectorTools::interpolate(dg->dof_handler,initial_condition,dg->solution);
+	// Create ODE solver using the factory and providing the DG object
+	std::shared_ptr<PHiLiP::ODE::ODESolver<dim, double>> ode_solver = PHiLiP::ODE::ODESolverFactory<dim, double>::create_ODESolver(dg);
+
+	double finalTime = 1.5;
+
+	//double dt = all_parameters->ode_solver_param.initial_time_step;
+	ode_solver->advance_solution_time(finalTime);
+
+	return 0; //need to change
+}
+
+int main (int argc, char * argv[])
+{
+	//parse parameters first
+	feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan
+	dealii::deallog.depth_console(99);
+		dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+		const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+		const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+		dealii::ConditionalOStream pcout(std::cout, mpi_rank==0);
+		pcout << "Starting program with " << n_mpi << " processors..." << std::endl;
+		if ((PHILIP_DIM==1) && !(n_mpi==1)) {
+			std::cout << "********************************************************" << std::endl;
+			std::cout << "Can't use mpirun -np X, where X>1, for 1D." << std::endl
+					  << "Currently using " << n_mpi << " processors." << std::endl
+					  << "Aborting..." << std::endl;
+			std::cout << "********************************************************" << std::endl;
+			std::abort();
+		}
+	int test_error = 1;
+	try
+	{
+        // Declare possible inputs
+        dealii::ParameterHandler parameter_handler;
+        PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler);
+        PHiLiP::Parameters::parse_command_line (argc, argv, parameter_handler);
+
+        // Read inputs from parameter file and set those values in AllParameters object
+        PHiLiP::Parameters::AllParameters all_parameters;
+        std::cout << "Reading input..." << std::endl;
+        all_parameters.parse_parameters (parameter_handler);
+
+        AssertDimension(all_parameters.dimension, PHILIP_DIM);
+
+        std::cout << "Starting program..." << std::endl;
+
+		using namespace PHiLiP;
+		//const Parameters::AllParameters parameters_input;
+		AdvectionPeriodic<PHILIP_DIM, 1> advection_test(&all_parameters);
+		int i = advection_test.run_test();
+		return i;
+	}
+	catch (std::exception &exc)
+	{
+		std::cerr << std::endl << std::endl
+				  << "----------------------------------------------------"
+				  << std::endl
+				  << "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
+	              << "Unknown exception!" << std::endl
+	              << "Aborting!" << std::endl
+	              << "----------------------------------------------------"
+	              << std::endl;
+	    return 1;
+	}
+	std::cout << "End of program." << std::endl;
+	return test_error;
+}
+
+
+
+
+
diff --git a/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm b/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm
index bad0610c8..490ea3952 100644
--- a/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm
+++ b/tests/burgers_inviscid_implicit/1d_burgers_inviscid_implicit.prm
@@ -5,7 +5,7 @@ set dimension = 1
 
 # The PDE we want to solve. Choices are
 # <diffusion|diffusion|convection_diffusion>.
-set pde_type  = burgers_inviscid
+set pde_type  = burgers_inviscid 
 
 subsection ODE solver
 
diff --git a/tests/burgers_inviscid_implicit/CMakeLists.txt b/tests/burgers_inviscid_implicit/CMakeLists.txt
index 5a42d1f51..239c9aa19 100644
--- a/tests/burgers_inviscid_implicit/CMakeLists.txt
+++ b/tests/burgers_inviscid_implicit/CMakeLists.txt
@@ -11,7 +11,7 @@ add_test(
 #add_test(
 #  NAME 2D_BURGERS_INVISCID_IMPLICIT_MANUFACTURED_SOLUTION
 #  COMMAND ${EXECUTABLE_OUTPUT_PATH}/PHiLiP_2D -i ${CMAKE_CURRENT_BINARY_DIR}/2d_burgers_inviscid_implicit.prm
-#  WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+#   WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
 #)
 
 #configure_file(3d_burgers_inviscid_implicit.prm 3d_burgers_inviscid_implicit.prm COPYONLY)
diff --git a/tests/burgers_split_stability/1D_burgers_stability.prm b/tests/burgers_split_stability/1D_burgers_stability.prm
new file mode 100644
index 000000000..18558e66d
--- /dev/null
+++ b/tests/burgers_split_stability/1D_burgers_stability.prm
@@ -0,0 +1,30 @@
+# -------------------
+# Number of dimensions
+set dimension = 1
+
+set use_weak_form = false
+
+set use_collocated_nodes = true
+
+set use_split_form = true
+
+set use_periodic_bc = true
+
+# The PDE we want to solve
+set pde_type = burgers_inviscid
+
+set conv_num_flux = split_form
+
+subsection ODE solver
+
+  set ode_output = verbose
+  
+  set nonlinear_max_iterations = 500
+
+  set print_iteration_modulo = 100
+
+  set ode_solver_type = explicit
+
+  set initial_time_step = 0.0001
+
+end
diff --git a/tests/burgers_split_stability/CMakeLists.txt b/tests/burgers_split_stability/CMakeLists.txt
new file mode 100644
index 000000000..9d7b75365
--- /dev/null
+++ b/tests/burgers_split_stability/CMakeLists.txt
@@ -0,0 +1,52 @@
+set(TEST_SRC
+    burgers_stability.cpp
+    )
+
+foreach(dim RANGE 1 1)
+
+    
+    # Output executable
+    string(CONCAT TEST_TARGET ${dim}D_burgers_energy_stability)
+    message("Adding executable " ${TEST_TARGET} " with files " ${TEST_SRC} "\n")
+    add_executable(${TEST_TARGET} ${TEST_SRC})
+    # Replace occurences of PHILIP_DIM with 1, 2, or 3 in the code
+    target_compile_definitions(${TEST_TARGET} PRIVATE PHILIP_DIM=${dim})
+
+    # Compile this executable when 'make unit_tests'
+    add_dependencies(unit_tests ${TEST_TARGET})
+    add_dependencies(${dim}D ${TEST_TARGET})
+
+    # Library dependency
+    set(ParameterLib ParametersLibrary)
+    string(CONCAT PhysicsLib Physics_${dim}D)
+    string(CONCAT NumericalFluxLib NumericalFlux_${dim}D)
+    string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) 
+    string(CONCAT ODESolverLib ODESolver_${dim}D)
+
+    target_link_libraries(${TEST_TARGET} ${ParameterLib})
+    target_link_libraries(${TEST_TARGET} ${PhysicsLib})
+    target_link_libraries(${TEST_TARGET} ${NumericalFluxLib})
+    target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib})
+    target_link_libraries(${TEST_TARGET} ${ODESolverLib})
+
+    target_link_libraries(${TEST_TARGET} m mpi mpi_cxx)
+
+    # Setup target with deal.II
+    DEAL_II_SETUP_TARGET(${TEST_TARGET})
+    
+    configure_file(1D_burgers_stability.prm 1D_burgers_stability.prm COPYONLY)
+    add_test(
+      NAME ${TEST_TARGET}
+      COMMAND ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} -i ${CMAKE_CURRENT_BINARY_DIR}/1D_burgers_stability.prm
+      WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+    )
+
+    unset(dim)
+    unset(TEST_TARGET)
+    unset(PhysicsLib)
+    unset(NumericalFluxLib)
+    unset(ParameterLib)
+    unset(DiscontinuousGalerkinLib)
+    unset(ODESolverLib)
+
+endforeach()
diff --git a/tests/burgers_split_stability/burgers_stability.cpp b/tests/burgers_split_stability/burgers_stability.cpp
new file mode 100644
index 000000000..ad98dae09
--- /dev/null
+++ b/tests/burgers_split_stability/burgers_stability.cpp
@@ -0,0 +1,212 @@
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/function.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/base/numbers.h>
+#include <deal.II/base/function_parser.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+
+#include "parameters/all_parameters.h"
+#include "parameters/parameters.h"
+#include "numerical_flux/numerical_flux.h"
+#include "physics/physics_factory.h"
+#include "physics/physics.h"
+#include "dg/dg.h"
+#include "ode_solver/ode_solver.h"
+
+#include<fenv.h>
+
+//using PDEType  = PHiLiP::Parameters::AllParameters::PartialDifferentialEquation;
+//using ConvType = PHiLiP::Parameters::AllParameters::ConvectiveNumericalFlux;
+//using DissType = PHiLiP::Parameters::AllParameters::DissipativeNumericalFlux;
+//
+//
+//const double TOLERANCE = 1E-12;
+
+
+template <int dim, int nstate>
+class BurgersEnergyStability
+{
+public:
+	BurgersEnergyStability() = delete;
+	BurgersEnergyStability(const PHiLiP::Parameters::AllParameters *const parameters_input);
+	int run_test();
+
+
+
+
+private:
+	double compute_energy(std::shared_ptr < PHiLiP::DGBase<dim, double> > &dg);
+    const PHiLiP::Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters
+};
+
+template <int dim, int nstate>
+BurgersEnergyStability<dim, nstate>::BurgersEnergyStability(const PHiLiP::Parameters::AllParameters *const parameters_input)
+:
+all_parameters(parameters_input)
+{}
+
+template<int dim, int nstate>
+double BurgersEnergyStability<dim, nstate>::compute_energy(std::shared_ptr < PHiLiP::DGBase<dim, double> > &dg)
+{
+	double energy = 0.0;
+	for (unsigned int i = 0; i < dg->solution.size(); ++i)
+	{
+		energy += 1./(dg->global_inverse_mass_matrix(i,i)) * dg->solution(i) * dg->solution(i);
+	}
+	return energy;
+}
+
+template <int dim, int nstate>
+int BurgersEnergyStability<dim, nstate>::run_test()
+{
+//	dealii::Triangulation<dim> grid(
+//	                typename dealii::Triangulation<dim>::MeshSmoothing(
+//	                    dealii::Triangulation<dim>::smoothing_on_refinement |
+//	                    dealii::Triangulation<dim>::smoothing_on_coarsening));
+	dealii::Triangulation<dim> grid;
+
+	double left = 0.0;
+	double right = 2.0;
+	const bool colorize = true;
+	int n_refinements = 5;
+	unsigned int poly_degree = 7;
+	dealii::GridGenerator::hyper_cube(grid, left, right, colorize);
+	grid.refine_global(n_refinements);
+	std::cout << "Grid generated and refined" << std::endl;
+	std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree);
+	std::cout << "dg created" <<std::endl;
+	dg->set_triangulation(&grid);
+	dg->allocate_system ();
+
+	std::cout << "Implement initial conditions" << std::endl;
+	dealii::FunctionParser<1> initial_condition;
+	std::string variables = "x";
+	std::map<std::string,double> constants;
+	constants["pi"] = dealii::numbers::PI;
+	std::string expression = "sin(pi*(x)) + 0.01";
+	initial_condition.initialize(variables,
+	                             expression,
+	                             constants);
+	dealii::VectorTools::interpolate(dg->dof_handler,initial_condition,dg->solution);
+	// Create ODE solver using the factory and providing the DG object
+	std::shared_ptr<PHiLiP::ODE::ODESolver<dim, double>> ode_solver = PHiLiP::ODE::ODESolverFactory<dim, double>::create_ODESolver(dg);
+
+	double finalTime = 3.;
+
+	double dt = all_parameters->ode_solver_param.initial_time_step;
+	//(void) dt;
+
+	//need to call ode_solver before calculating energy because mass matrix isn't allocated yet.
+
+	ode_solver->advance_solution_time(0.000001);
+	double initial_energy = compute_energy(dg);
+
+	//currently the only way to calculate energy at each time-step is to advance solution by dt instead of finaltime
+	//this causes some issues with outputs (only one file is output, which is overwritten at each time step)
+	//also the ode solver output doesn't make sense (says "iteration 1 out of 1")
+	//but it works. I'll keep it for now and need to modify the output functions later to account for this.
+	std::ofstream myfile ("energy_plot.gpl" , std::ios::trunc);
+
+	for (int i = 0; i < std::ceil(finalTime/dt); ++ i)
+	{
+		ode_solver->advance_solution_time(dt);
+		double current_energy = compute_energy(dg);
+		std::cout << "Energy at time " << i * dt << " is " << current_energy << std::endl;
+		myfile << i * dt << " " << current_energy << std::endl;
+		if (current_energy - initial_energy >= 0.001)
+		{
+			return 1;
+			break;
+		}
+	}
+	myfile.close();
+
+
+	//ode_solver->advance_solution_time(finalTime);
+
+	return 0; //need to change
+}
+
+int main (int argc, char * argv[])
+{
+	//parse parameters first
+	feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan
+
+	dealii::deallog.depth_console(99);
+	dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+	const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+	const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+	dealii::ConditionalOStream pcout(std::cout, mpi_rank==0);
+	pcout << "Starting program with " << n_mpi << " processors..." << std::endl;
+	if ((PHILIP_DIM==1) && !(n_mpi==1)) {
+		std::cout << "********************************************************" << std::endl;
+		std::cout << "Can't use mpirun -np X, where X>1, for 1D." << std::endl
+				  << "Currently using " << n_mpi << " processors." << std::endl
+				  << "Aborting..." << std::endl;
+		std::cout << "********************************************************" << std::endl;
+		std::abort();
+	}
+
+	int test_error = 1;
+	try
+	{
+        // Declare possible inputs
+        dealii::ParameterHandler parameter_handler;
+        PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler);
+        PHiLiP::Parameters::parse_command_line (argc, argv, parameter_handler);
+
+        // Read inputs from parameter file and set those values in AllParameters object
+        PHiLiP::Parameters::AllParameters all_parameters;
+        std::cout << "Reading input..." << std::endl;
+        all_parameters.parse_parameters (parameter_handler);
+
+        AssertDimension(all_parameters.dimension, PHILIP_DIM);
+
+        std::cout << "Starting program..." << std::endl;
+
+		//using namespace PHiLiP;
+		//const Parameters::AllParameters parameters_input;
+		BurgersEnergyStability<PHILIP_DIM, PHILIP_DIM> burgers_test(&all_parameters);
+		std::cout << "Running test" <<std::endl;
+		int i = burgers_test.run_test();
+		return i;
+	}
+	catch (std::exception &exc)
+	{
+		std::cerr << std::endl << std::endl
+				  << "----------------------------------------------------"
+				  << std::endl
+				  << "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
+	              << "Unknown exception!" << std::endl
+	              << "Aborting!" << std::endl
+	              << "----------------------------------------------------"
+	              << std::endl;
+	    return 1;
+	}
+	std::cout << "End of program." << std::endl;
+	return test_error;
+}
+
+
+
+
+
diff --git a/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm b/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm
index 124a69394..71a33fabb 100644
--- a/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm
+++ b/tests/euler_integration/1d_euler_laxfriedrichs_manufactured.prm
@@ -7,6 +7,7 @@ set pde_type  = euler
 
 set conv_num_flux  = lax_friedrichs
 
+
 subsection ODE solver
 
   set ode_output                          = verbose
diff --git a/tests/euler_integration/2d_euler_laxfriedrichs_manufactured.prm b/tests/euler_integration/2d_euler_laxfriedrichs_manufactured.prm
index 12e3f1dc3..f8caab53e 100644
--- a/tests/euler_integration/2d_euler_laxfriedrichs_manufactured.prm
+++ b/tests/euler_integration/2d_euler_laxfriedrichs_manufactured.prm
@@ -7,6 +7,7 @@ set pde_type  = euler
 
 set conv_num_flux  = lax_friedrichs
 
+
 subsection ODE solver
   # Maximum nonlinear solver iterations
   set nonlinear_max_iterations            = 100
diff --git a/tests/euler_integration/2d_euler_vortex.prm b/tests/euler_integration/2d_euler_vortex.prm
index 460a71178..697e78782 100644
--- a/tests/euler_integration/2d_euler_vortex.prm
+++ b/tests/euler_integration/2d_euler_vortex.prm
@@ -10,6 +10,7 @@ set pde_type  = euler
 
 set conv_num_flux = roe
 
+
 subsection euler
   set reference_length = 1.0
   set mach_infinity = 0.5
diff --git a/tests/euler_split_inviscid_taylor_green_vortex/3D_euler_split_inviscid_taylor_green_vortex.prm b/tests/euler_split_inviscid_taylor_green_vortex/3D_euler_split_inviscid_taylor_green_vortex.prm
new file mode 100644
index 000000000..a3946e579
--- /dev/null
+++ b/tests/euler_split_inviscid_taylor_green_vortex/3D_euler_split_inviscid_taylor_green_vortex.prm
@@ -0,0 +1,30 @@
+# -------------------
+# Number of dimensions
+set dimension = 3
+
+set use_weak_form = false
+
+set use_collocated_nodes = true
+
+set use_split_form = false
+
+set use_periodic_bc = true
+
+# The PDE we want to solve
+set pde_type = euler
+
+set conv_num_flux = split_form
+
+subsection ODE solver
+
+  set ode_output = verbose
+  
+  set nonlinear_max_iterations = 500
+
+  set print_iteration_modulo = 10
+
+  set ode_solver_type = explicit
+
+  set initial_time_step = 0.001
+
+end
diff --git a/tests/euler_split_inviscid_taylor_green_vortex/CMakeLists.txt b/tests/euler_split_inviscid_taylor_green_vortex/CMakeLists.txt
new file mode 100644
index 000000000..b5754c2f9
--- /dev/null
+++ b/tests/euler_split_inviscid_taylor_green_vortex/CMakeLists.txt
@@ -0,0 +1,52 @@
+set(TEST_SRC
+    euler_split_inviscid_taylor_green_vortex.cpp
+    )
+
+foreach(dim RANGE 3 3)
+
+    
+    # Output executable
+    string(CONCAT TEST_TARGET ${dim}D_euler_split_inviscid_taylor_green_vortex)
+    message("Adding executable " ${TEST_TARGET} " with files " ${TEST_SRC} "\n")
+    add_executable(${TEST_TARGET} ${TEST_SRC})
+    # Replace occurences of PHILIP_DIM with 1, 2, or 3 in the code
+    target_compile_definitions(${TEST_TARGET} PRIVATE PHILIP_DIM=${dim})
+
+    # Compile this executable when 'make unit_tests'
+    add_dependencies(unit_tests ${TEST_TARGET})
+    add_dependencies(${dim}D ${TEST_TARGET})
+
+    # Library dependency
+    set(ParameterLib ParametersLibrary)
+    string(CONCAT PhysicsLib Physics_${dim}D)
+    string(CONCAT NumericalFluxLib NumericalFlux_${dim}D)
+    string(CONCAT DiscontinuousGalerkinLib DiscontinuousGalerkin_${dim}D) 
+    string(CONCAT ODESolverLib ODESolver_${dim}D)
+
+    target_link_libraries(${TEST_TARGET} ${ParameterLib})
+    target_link_libraries(${TEST_TARGET} ${PhysicsLib})
+    target_link_libraries(${TEST_TARGET} ${NumericalFluxLib})
+    target_link_libraries(${TEST_TARGET} ${DiscontinuousGalerkinLib})
+    target_link_libraries(${TEST_TARGET} ${ODESolverLib})
+
+    #target_link_libraries(${TEST_TARGET} m mpi mpi_cxx)
+
+    # Setup target with deal.II
+    DEAL_II_SETUP_TARGET(${TEST_TARGET})
+    
+    configure_file(3D_euler_split_inviscid_taylor_green_vortex.prm 3D_euler_split_inviscid_taylor_green_vortex.prm COPYONLY)
+    add_test(
+      NAME ${TEST_TARGET}
+      COMMAND ${EXECUTABLE_OUTPUT_PATH}/${TEST_TARGET} -i ${CMAKE_CURRENT_BINARY_DIR}/3D_euler_split_inviscid_taylor_green_vortex.prm
+      WORKING_DIRECTORY ${TEST_OUTPUT_DIR}
+    )
+
+    unset(dim)
+    unset(TEST_TARGET)
+    unset(PhysicsLib)
+    unset(NumericalFluxLib)
+    unset(ParameterLib)
+    unset(DiscontinuousGalerkinLib)
+    unset(ODESolverLib)
+
+endforeach()
diff --git a/tests/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp b/tests/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp
new file mode 100644
index 000000000..9d104faee
--- /dev/null
+++ b/tests/euler_split_inviscid_taylor_green_vortex/euler_split_inviscid_taylor_green_vortex.cpp
@@ -0,0 +1,289 @@
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/function.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/solution_transfer.h>
+#include <deal.II/base/numbers.h>
+#include <deal.II/base/function_parser.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+
+#include "parameters/all_parameters.h"
+#include "parameters/parameters.h"
+#include "numerical_flux/numerical_flux.h"
+#include "physics/physics_factory.h"
+#include "physics/physics.h"
+#include "dg/dg.h"
+#include "ode_solver/ode_solver.h"
+
+#include<fenv.h>
+
+//using PDEType  = PHiLiP::Parameters::AllParameters::PartialDifferentialEquation;
+//using ConvType = PHiLiP::Parameters::AllParameters::ConvectiveNumericalFlux;
+//using DissType = PHiLiP::Parameters::AllParameters::DissipativeNumericalFlux;
+//
+//
+//const double TOLERANCE = 1E-12;
+
+
+template <int dim, int nstate>
+class EulerTaylorGreen
+{
+public:
+	EulerTaylorGreen() = delete;
+	EulerTaylorGreen(const PHiLiP::Parameters::AllParameters *const parameters_input);
+	int run_test();
+
+
+
+
+private:
+	double compute_kinetic_energy(std::shared_ptr < PHiLiP::DGBase<dim, double> > &dg, unsigned int poly_degree);
+	double compute_quadrature_kinetic_energy(std::array<double,nstate> soln_at_q);
+    const PHiLiP::Parameters::AllParameters *const all_parameters; ///< Pointer to all parameters
+    const MPI_Comm mpi_communicator;
+    dealii::ConditionalOStream pcout;
+};
+
+template <int dim, int nstate>
+EulerTaylorGreen<dim, nstate>::EulerTaylorGreen(const PHiLiP::Parameters::AllParameters *const parameters_input)
+:
+all_parameters(parameters_input)
+, mpi_communicator(MPI_COMM_WORLD)
+, pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0)
+{}
+
+template <int dim, int nstate>
+double EulerTaylorGreen<dim,nstate>::compute_quadrature_kinetic_energy(std::array<double,nstate> soln_at_q)
+{
+	const double density = soln_at_q[0];
+
+	const double quad_kin_energ =  0.5*(soln_at_q[1]*soln_at_q[1] +
+										soln_at_q[2]*soln_at_q[2] +
+										soln_at_q[3]*soln_at_q[3] )/density;
+	return quad_kin_energ;
+}
+
+template<int dim, int nstate>
+double EulerTaylorGreen<dim, nstate>::compute_kinetic_energy(std::shared_ptr < PHiLiP::DGBase<dim, double> > &dg, unsigned int poly_degree)
+{
+	// Overintegrate the error to make sure there is not integration error in the error estimate
+	int overintegrate = 10 ;//10;
+	dealii::QGauss<dim> quad_extra(dg->max_degree+1+overintegrate);
+	            dealii::FEValues<dim,dim> fe_values_extra(dealii::MappingQ<dim>(dg->max_degree+overintegrate), dg->fe_collection[poly_degree], quad_extra,
+	            dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
+//	dealii::QGauss<dim> quad_extra(dg->fe_system.tensor_degree()+overintegrate);
+//	dealii::FEValues<dim,dim> fe_values_extra(dg->mapping, dg->fe_system, quad_extra,
+//							dealii::update_values | dealii::update_JxW_values | dealii::update_quadrature_points);
+	const unsigned int n_quad_pts = fe_values_extra.n_quadrature_points;
+	std::array<double,nstate> soln_at_q;
+
+	double total_kinetic_energy = 0;
+
+	// Integrate solution error and output error
+	typename dealii::DoFHandler<dim>::active_cell_iterator
+	cell = dg->dof_handler.begin_active(),
+	endc = dg->dof_handler.end();
+
+	std::vector<dealii::types::global_dof_index> dofs_indices (fe_values_extra.dofs_per_cell);
+
+	            //const double gam = euler_physics_double.gam;
+	            //const double mach_inf = euler_physics_double.mach_inf;
+	            //const double tot_temperature_inf = 1.0;
+	            //const double tot_pressure_inf = 1.0;
+	            //// Assuming a tank at rest, velocity = 0, therefore, static pressure and temperature are same as total
+	            //const double density_inf = gam*tot_pressure_inf/tot_temperature_inf * mach_inf * mach_inf;
+	            //const double entropy_inf = tot_pressure_inf*pow(density_inf,-gam);
+	//const double entropy_inf = euler_physics_double.entropy_inf;
+
+	for (; cell!=endc; ++cell) {
+
+		fe_values_extra.reinit (cell);
+		//std::cout << "sitting on cell " << cell->index() << std::endl;
+	    cell->get_dof_indices (dofs_indices);
+        for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
+
+        	std::fill(soln_at_q.begin(), soln_at_q.end(), 0);
+	        for (unsigned int idof=0; idof<fe_values_extra.dofs_per_cell; ++idof) {
+	        	const unsigned int istate = fe_values_extra.get_fe().system_to_component_index(idof).first;
+	            soln_at_q[istate] += dg->solution[dofs_indices[idof]] * fe_values_extra.shape_value_component(idof, iquad, istate);
+	        }
+	        const double quadrature_kinetic_energy = compute_quadrature_kinetic_energy(soln_at_q);
+
+	        total_kinetic_energy += quadrature_kinetic_energy * fe_values_extra.JxW(iquad);
+        }
+	}
+	return total_kinetic_energy;
+}
+
+template <int dim, int nstate>
+int EulerTaylorGreen<dim, nstate>::run_test()
+{
+	//dealii::Triangulation<dim> grid;
+	dealii::parallel::distributed::Triangulation<dim> grid(this->mpi_communicator);
+
+	double left = 0.0;
+	double right = 2 * dealii::numbers::PI;
+	const bool colorize = true;
+	int n_refinements = 2;
+	unsigned int poly_degree = 1;
+	dealii::GridGenerator::hyper_cube(grid, left, right, colorize);
+
+	std::vector<dealii::GridTools::PeriodicFacePair<typename dealii::Triangulation<PHILIP_DIM>::cell_iterator> > matched_pairs;
+	dealii::GridTools::collect_periodic_faces(grid,0,1,0,matched_pairs);
+	dealii::GridTools::collect_periodic_faces(grid,2,3,1,matched_pairs);
+	dealii::GridTools::collect_periodic_faces(grid,4,5,2,matched_pairs);
+	grid.add_periodicity(matched_pairs);
+
+	grid.refine_global(n_refinements);
+
+	std::shared_ptr < PHiLiP::DGBase<dim, double> > dg = PHiLiP::DGFactory<dim,double>::create_discontinuous_galerkin(all_parameters, poly_degree);
+	dg->set_triangulation(&grid);
+	dg->allocate_system ();
+
+	std::cout << "Implement initial conditions" << std::endl;
+	dealii::FunctionParser<PHILIP_DIM> initial_condition(5);
+	std::string variables = "x,y,z";
+	std::map<std::string,double> constants;
+	constants["pi"] = dealii::numbers::PI;
+	std::vector<std::string> expressions(5);
+	expressions[0] = "1";
+	expressions[1] = "sin(x)*cos(y)*cos(z)";
+	expressions[2] = "-cos(x)*sin(y)*cos(z)";
+	expressions[3] = "0";
+	expressions[4] = "250.0/1.4 + 2.5/16.0 * (cos(2.0*x)*cos(2.0*z) + 2.0*cos(2.0*y) + 2.0*cos(2.0*x) + cos(2.0*y)*cos(2.0*z)) + 0.5 * pow(cos(z),2.0) * (pow(cos(x),2.0) * pow(sin(y),2.0) +pow(sin(x),2.0) * pow(cos(y),2.0))";
+	initial_condition.initialize(variables,
+	                             expressions,
+	                             constants);
+//	dealii::Point<3> point (0.0,1.0,1.0);
+//	dealii::Vector<double> result(5);
+//	initial_condition.vector_value(point,result);
+//	dealii::deallog << result;
+
+	std::cout << "initial condition successfully implemented" << std::endl;
+	dealii::VectorTools::interpolate(dg->dof_handler,initial_condition,dg->solution);
+	std::cout << "initial condition interpolated to DG solution" << std::endl;
+	// Create ODE solver using the factory and providing the DG object
+
+	std::cout << "creating ODE solver" << std::endl;
+	std::shared_ptr<PHiLiP::ODE::ODESolver<dim, double>> ode_solver = PHiLiP::ODE::ODESolverFactory<dim, double>::create_ODESolver(dg);
+	std::cout << "ODE solver successfully created" << std::endl;
+	double finalTime = 14.;
+	double dt = all_parameters->ode_solver_param.initial_time_step;
+
+	//double dt = all_parameters->ode_solver_param.initial_time_step;
+	//need to call ode_solver before calculating energy because mass matrix isn't allocated yet.
+	//ode_solver->advance_solution_time(0.000001);
+	//double initial_energy = compute_energy(dg);
+
+	std::cout << "preparing to advance solution in time" << std::endl;
+//	//(void) finalTime;
+//	std::cout << "kinetic energy is" << compute_kinetic_energy(dg) <<std::endl;
+//	ode_solver->advance_solution_time(finalTime);
+//	//std::cout << "kinetic energy is" << compute_kinetic_energy(dg) <<std::endl;
+//	std::cout << "kinetic energy is" << compute_kinetic_energy(dg) <<std::endl;
+
+
+	//currently the only way to calculate energy at each time-step is to advance solution by dt instead of finaltime
+	//this causes some issues with outputs (only one file is output, which is overwritten at each time step)
+	//also the ode solver output doesn't make sense (says "iteration 1 out of 1")
+	//but it works. I'll keep it for now and need to modify the output functions later to account for this.
+
+	std::ofstream myfile ("kinetic_energy_plot.gpl" , std::ios::trunc);
+
+	for (int i = 0; i < std::ceil(finalTime/dt); ++ i)
+	{
+		ode_solver->advance_solution_time(dt);
+		double current_energy = compute_kinetic_energy(dg,poly_degree);
+		std::cout << "Energy at time " << i * dt << " is " << current_energy << std::endl;
+		myfile << i * dt << " " << current_energy << std::endl;
+
+	}
+
+	myfile.close();
+
+//	ode_solver->advance_solution_time(finalTime);
+//	(void) dt;
+	return 0; //need to change
+}
+
+int main (int argc, char * argv[])
+{
+	//parse parameters first
+	feenableexcept(FE_INVALID | FE_OVERFLOW); // catch nan
+	dealii::deallog.depth_console(99);
+		dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+		const int n_mpi = dealii::Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+		const int mpi_rank = dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+		dealii::ConditionalOStream pcout(std::cout, mpi_rank==0);
+		pcout << "Starting program with " << n_mpi << " processors..." << std::endl;
+		if ((PHILIP_DIM==1) && !(n_mpi==1)) {
+			std::cout << "********************************************************" << std::endl;
+			std::cout << "Can't use mpirun -np X, where X>1, for 1D." << std::endl
+					  << "Currently using " << n_mpi << " processors." << std::endl
+					  << "Aborting..." << std::endl;
+			std::cout << "********************************************************" << std::endl;
+			std::abort();
+		}
+	int test_error = 1;
+	try
+	{
+        // Declare possible inputs
+        dealii::ParameterHandler parameter_handler;
+        PHiLiP::Parameters::AllParameters::declare_parameters (parameter_handler);
+        PHiLiP::Parameters::parse_command_line (argc, argv, parameter_handler);
+
+        // Read inputs from parameter file and set those values in AllParameters object
+        PHiLiP::Parameters::AllParameters all_parameters;
+        std::cout << "Reading input..." << std::endl;
+        all_parameters.parse_parameters (parameter_handler);
+
+        AssertDimension(all_parameters.dimension, PHILIP_DIM);
+
+        std::cout << "Starting program..." << std::endl;
+
+		using namespace PHiLiP;
+		//const Parameters::AllParameters parameters_input;
+		EulerTaylorGreen<PHILIP_DIM, PHILIP_DIM+2> euler_test(&all_parameters);
+		int i = euler_test.run_test();
+		return i;
+	}
+	catch (std::exception &exc)
+	{
+		std::cerr << std::endl << std::endl
+				  << "----------------------------------------------------"
+				  << std::endl
+				  << "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
+	              << "Unknown exception!" << std::endl
+	              << "Aborting!" << std::endl
+	              << "----------------------------------------------------"
+	              << std::endl;
+	    return 1;
+	}
+	std::cout << "End of program." << std::endl;
+	return test_error;
+}
+
+
+
+
+