From 593089b2a5de8f540c4c68efea42d2b6c6742403 Mon Sep 17 00:00:00 2001
From: dougshidong <doug.shidong@gmail.com>
Date: Fri, 17 May 2019 21:13:33 -0400
Subject: [PATCH] Major commit taking out numerical flux

Physics and numerical fluxes are now computed outside of the DG discretization.
Partial implementation of vector-valued problems by adding the std::array structures

Diffusive terms require a function that will evaluate the solution numerical flux and the viscous numerical flux.
Current implementation of SIPG is working to recover optimal convergence orders of the solution but not of the output functional

Suboptimal output error is perceived for diffusion (maybe, I get 2p) and convection-diffusion (definitely)

TODO: Find the adjoint consistent boundary condition for the diffusion problem and convection-diffusion problem
TODO: Find out why the linear solver is soooo much slower for the diffusion problems

      Start  1: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION
 1/12 Test  #1: 1D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ..............   Passed    0.53 sec
      Start  2: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION
 2/12 Test  #2: 2D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ..............   Passed    5.84 sec
      Start  3: 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION
 3/12 Test  #3: 3D_ADVECTION_IMPLICIT_MANUFACTURED_SOLUTION ..............   Passed   47.98 sec
      Start  4: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION
 4/12 Test  #4: 1D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ..............   Passed    0.64 sec
      Start  5: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION
 5/12 Test  #5: 2D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ..............   Passed    3.23 sec
      Start  6: 3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION
 6/12 Test  #6: 3D_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ..............   Passed  237.84 sec
      Start  7: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION
 7/12 Test  #7: 1D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ...   Passed    0.75 sec
      Start  8: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION
 8/12 Test  #8: 2D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ...   Passed    3.45 sec
      Start  9: 3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION
 9/12 Test  #9: 3D_CONVECTION_DIFFUSION_IMPLICIT_MANUFACTURED_SOLUTION ...   Passed  266.84 sec
      Start 10: 1D_numerical_flux_conservation
10/12 Test #10: 1D_numerical_flux_conservation ...........................   Passed    0.27 sec
      Start 11: 2D_numerical_flux_conservation
11/12 Test #11: 2D_numerical_flux_conservation ...........................   Passed    0.25 sec
      Start 12: 3D_numerical_flux_conservation
12/12 Test #12: 3D_numerical_flux_conservation ...........................   Passed    0.24 sec

100% tests passed, 0 tests failed out of 12

Total Test time (real) = 567.87 sec
---
 TODO                                          |  10 +-
 input_file.prm                                |  43 +++
 run_this_input.py                             |  28 ++
 src/assemble_implicit.cpp                     | 219 ++++++------
 src/dg.cpp                                    |  41 +++
 src/dg.h                                      |   4 +
 src/grid_study.cpp                            |   4 +-
 src/linear_solver.cpp                         |  13 +-
 src/numerical_flux/CMakeLists.txt             |   1 +
 src/numerical_flux/numerical_flux.cpp         |  30 +-
 src/numerical_flux/numerical_flux.h           |  51 +--
 src/numerical_flux/viscous_numerical_flux.cpp | 157 +++++++++
 src/numerical_flux/viscous_numerical_flux.h   |  77 +++++
 src/ode_solver.cpp                            |  22 +-
 src/parameters.cpp                            |   2 +-
 src/physics/CMakeLists.txt                    |   6 +-
 src/physics/physics.cpp                       | 312 ++++++++++++++++++
 src/physics/physics.h                         |  33 +-
 .../1d_convection_diffusion_implicit.prm      |  13 +-
 .../2d_convection_diffusion_implicit.prm      |   6 +-
 .../3d_convection_diffusion_implicit.prm      |   4 +-
 .../2d_diffusion_implicit.prm                 |   4 +-
 .../numerical_flux_conservation.cpp           | 209 ++++++++----
 23 files changed, 1000 insertions(+), 289 deletions(-)
 create mode 100644 input_file.prm
 create mode 100755 run_this_input.py
 create mode 100644 src/numerical_flux/viscous_numerical_flux.cpp
 create mode 100644 src/numerical_flux/viscous_numerical_flux.h

diff --git a/TODO b/TODO
index fd51957e9..b78031f38 100644
--- a/TODO
+++ b/TODO
@@ -2,15 +2,11 @@ Look into member variable naming for fe_values
     This one just bit me.
     In face_boundary I used the member fe_values instead of passed fe_values_face by mistake.
 
-Take physics out of the DiscontinuousGalerkin class.
+Find the adjoint consistent boundary condition for the diffusion problem and convection-diffusion problem
 
-Move boundary conditions to Physics class.
-
-Compute auxiliary variable $ \sigma = \nabla u $ independently.
+Find out why the linear solver is soooo much slower for the diffusion problems. Maybe just need pre-conditioner.
 
-Take numerical flux out of the DiscontinuousGalerkin class.
-
-Adjoint consistency, output integral 2p.
+Move boundary conditions to Physics class.
 
 Vector-valued problems
     - Figure out global solution ordering
diff --git a/input_file.prm b/input_file.prm
new file mode 100644
index 000000000..4d98bf789
--- /dev/null
+++ b/input_file.prm
@@ -0,0 +1,43 @@
+# Listing of Parameters
+# ---------------------
+# Number of dimensions
+set dimension = 1
+
+# The kind of solver for the linear system. Choices are
+# <advection|convection_diffusion>.
+set pde_type  = advection
+
+
+subsection ODE solver
+  # Maximum nonlinear solver iterations
+  set nonlinear_max_iterations            = 3
+
+  # Nonlinear solver residual tolerance
+  set nonlinear_steady_residual_tolerance = 1e-13
+
+  # Print every print_iteration_modulo iterations of the nonlinear solver
+  set print_iteration_modulo              = 1
+
+  # Explicit or implicit solverChoices are <explicit|implicit>.
+  set solver_type                         = implicit
+end
+
+
+subsection manufactured solution convergence study
+  # Last degree used for convergence study
+  set degree_end   = 4
+
+  # Starting degree for convergence study
+  set degree_start = 1
+
+  # Initial grid of size (initial_grid_size)^dim
+  set initial_grid_size = 2
+
+  # Number of grids in grid study
+  set number_of_grids   = 5
+
+  # Multiplier on grid size. nth-grid will be of size
+  # (initial_grid^grid_progression)^dim
+  set grid_progression  = 1.5
+end
+
diff --git a/run_this_input.py b/run_this_input.py
new file mode 100755
index 000000000..f46d9195f
--- /dev/null
+++ b/run_this_input.py
@@ -0,0 +1,28 @@
+#!/usr/bin/env python3
+import os
+import sys
+import re
+
+fname = 'input_file.prm'
+
+f = open(fname)
+
+regexp = re.compile(r'set dimension.*?([1-3.-]+)')
+
+dim = -1
+with open(fname) as f:
+    for line in f:
+        match = regexp.match(line)
+        if match:
+            dim = int(match.group(1))
+            break
+if dim == -1:
+    sys.exit("No valid 'set dimension = [1-3]' line found")
+
+
+
+bin_path = './bin/PHiLiP_'+str(dim)+'D'
+
+command = bin_path + ' -p ' + fname
+print("Running: " + command)
+os.system(command)
diff --git a/src/assemble_implicit.cpp b/src/assemble_implicit.cpp
index 2e938b85c..10d09b33b 100644
--- a/src/assemble_implicit.cpp
+++ b/src/assemble_implicit.cpp
@@ -26,7 +26,7 @@ namespace PHiLiP
     {
         using ADtype = Sacado::Fad::DFad<real>;
         using ADArray = std::array<ADtype,nstate>;
-        using ADArrayVector = std::array< Tensor<1,dim,ADtype>, nstate >;
+        using ADArrayTensor1 = std::array< Tensor<1,dim,ADtype>, nstate >;
 
         const unsigned int n_quad_pts      = fe_values_vol->n_quadrature_points;
         const unsigned int n_dofs_cell     = fe_values_vol->dofs_per_cell;
@@ -45,10 +45,10 @@ namespace PHiLiP
         std::vector<real> residual_derivatives(n_dofs_cell);
 
         std::vector< ADArray > soln_at_q(n_quad_pts);
-        std::vector< ADArrayVector > soln_grad_at_q(n_quad_pts); // Tensor initialize with zeros
+        std::vector< ADArrayTensor1 > soln_grad_at_q(n_quad_pts); // Tensor initialize with zeros
 
-        std::vector< ADArrayVector > conv_phys_flux_at_q(n_quad_pts);
-        std::vector< ADArrayVector > diss_phys_flux_at_q(n_quad_pts);
+        std::vector< ADArrayTensor1 > conv_phys_flux_at_q(n_quad_pts);
+        std::vector< ADArrayTensor1 > diss_phys_flux_at_q(n_quad_pts);
         std::vector< ADArray > source_at_q(n_quad_pts);
 
         for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
@@ -85,12 +85,12 @@ namespace PHiLiP
 
                 for (int istate=0; istate<nstate; istate++) {
                     // Convective
-                    rhs[istate] += conv_phys_flux_at_q[iquad][istate] * fe_values_vol->shape_grad(itest,iquad) * JxW[iquad];
+                    rhs[istate] += fe_values_vol->shape_grad(itest,iquad) * conv_phys_flux_at_q[iquad][istate] * JxW[iquad];
                     // Diffusive
                     // Note that for diffusion, the negative is defined in the physics
-                    //rhs[istate] += diss_phys_flux_at_q[iquad][istate] * fe_values_vol->shape_grad(itest,iquad) * JxW[iquad];
+                    rhs[istate] += fe_values_vol->shape_grad(itest,iquad) * diss_phys_flux_at_q[iquad][istate] * JxW[iquad];
                     // Source
-                    rhs[istate] += source_at_q[iquad][istate] * fe_values_vol->shape_value(itest,iquad) * JxW[iquad];
+                    rhs[istate] += fe_values_vol->shape_value(itest,iquad) * source_at_q[iquad][istate] * JxW[iquad];
                 }
             }
 
@@ -119,7 +119,7 @@ namespace PHiLiP
     {
         using ADtype = Sacado::Fad::DFad<real>;
         using ADArray = std::array<ADtype,nstate>;
-        using ADArrayVector = std::array< Tensor<1,dim,ADtype>, nstate >;
+        using ADArrayTensor1 = std::array< Tensor<1,dim,ADtype>, nstate >;
 
         const unsigned int n_dofs_cell = fe_values_boundary->dofs_per_cell;
         const unsigned int n_face_quad_pts = fe_values_boundary->n_quadrature_points;
@@ -136,10 +136,12 @@ namespace PHiLiP
         //boundary_function.value_list (fe_values_boundary->get_quadrature_points(), boundary_values, dummy);
 
         std::vector<ADArray> boundary_values(n_face_quad_pts);
+        std::vector<ADArrayTensor1> boundary_gradients(n_face_quad_pts);
         const std::vector< Point<dim,real> > quad_pts = fe_values_boundary->get_quadrature_points();
         for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
             const Point<dim, real> x_quad = quad_pts[iquad];
             pde_physics->manufactured_solution(x_quad, boundary_values[iquad]);
+            pde_physics->manufactured_gradient(x_quad, boundary_gradients[iquad]);
         }
 
         // AD variable
@@ -154,18 +156,22 @@ namespace PHiLiP
         std::vector<ADArray> soln_int(n_face_quad_pts);
         std::vector<ADArray> soln_ext(n_face_quad_pts);
 
-        std::vector<ADArrayVector> soln_grad_int(n_face_quad_pts);
-        std::vector<ADArrayVector> soln_grad_ext(n_face_quad_pts);
+        std::vector<ADArrayTensor1> soln_grad_int(n_face_quad_pts);
+        std::vector<ADArrayTensor1> soln_grad_ext(n_face_quad_pts);
 
-        std::vector<ADArrayVector> diss_phys_flux_int(n_face_quad_pts);
-        std::vector<ADArrayVector> diss_phys_flux_ext(n_face_quad_pts);
+        std::vector<ADArrayTensor1> diss_phys_flux_int(n_face_quad_pts);
+        std::vector<ADArrayTensor1> diss_phys_flux_ext(n_face_quad_pts);
 
 
-        std::vector<ADArray> convective_num_flux_dot_n(n_face_quad_pts);
+        std::vector<ADArray> conv_num_flux_dot_n(n_face_quad_pts);
+        std::vector<ADArray> diss_soln_num_flux(n_face_quad_pts); // u*
+        std::vector<ADArrayTensor1> diss_flux_jump_int(n_face_quad_pts); // u*-u_int
+        std::vector<ADArray> diss_auxi_num_flux_dot_n(n_face_quad_pts); // sigma*
 
         for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
 
             const Tensor<1,dim,ADtype> normal_int = normals[iquad];
+            const Tensor<1,dim,ADtype> normal_ext = -normal_int;
             ADArray characteristic_dot_n_at_q = pde_physics->convective_eigenvalues(boundary_values[iquad], normal_int);
 
             for (int istate=0; istate<nstate; istate++) {
@@ -180,18 +186,42 @@ namespace PHiLiP
 
                 const bool inflow = (characteristic_dot_n_at_q[istate] < 0.);
 
+                // For some reason, Nitsche boundary type is adjoint inconsistent at the inflow
+                // but we need to impose dirichlet Nitsche boundary type at outflows
                 if (inflow) { // Dirichlet boundary condition
                     soln_ext[iquad][istate] = boundary_values[iquad][istate];
-                    // For some reason, Nitsche boundary type is adjoint inconsistent
                     //soln_ext[iquad][istate] = -soln_int[iquad][istate]+2*boundary_values[iquad][istate];
+
+                    soln_grad_ext[iquad][istate] = soln_grad_int[iquad][istate];
                 } else { // Neumann boundary condition
-                    soln_ext[iquad][istate] = soln_int[iquad][istate];
+                    //soln_ext[iquad][istate] = soln_int[iquad][istate];
+                    //soln_ext[iquad][istate] = boundary_values[iquad][istate];
+                    soln_ext[iquad][istate] = -soln_int[iquad][istate]+2*boundary_values[iquad][istate];
+
+                    //soln_grad_ext[iquad][istate] = soln_grad_int[iquad][istate];
+                    soln_grad_ext[iquad][istate] = boundary_gradients[iquad][istate];
                 }
+
             }
-            // Evaluate physical convective flux and source term
+            // Evaluate physical convective flux, physical dissipative flux, and source term
             pde_physics->dissipative_flux (soln_int[iquad], soln_grad_int[iquad], diss_phys_flux_int[iquad]);
 
-            convective_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+            conv_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+
+            diss_soln_num_flux[iquad] = diss_num_flux->evaluate_solution_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+
+            ADArrayTensor1 diss_soln_jump_int;
+            for (int s=0; s<nstate; s++) {
+                diss_soln_jump_int[s] = (diss_soln_num_flux[iquad][s] - soln_int[iquad][s]) * normal_int;
+            }
+            pde_physics->dissipative_flux (soln_int[iquad], diss_soln_jump_int, diss_flux_jump_int[iquad]);
+
+            diss_auxi_num_flux_dot_n[iquad] = diss_num_flux->evaluate_auxiliary_flux(
+                soln_int[iquad], soln_ext[iquad],
+                soln_grad_int[iquad], soln_grad_ext[iquad],
+                normal_int, penalty);
+
+
         }
 
         // Applying convection boundary condition
@@ -205,11 +235,18 @@ namespace PHiLiP
             for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
 
                 for (int istate=0; istate<nstate; istate++) {
-                    rhs[istate] -= 
-                        convective_num_flux_dot_n[iquad][istate] *
-                        fe_values_boundary->shape_value(itest,iquad) *
-                        JxW[iquad];
+                    // Convection
+                    rhs[istate] -= fe_values_boundary->shape_value(itest,iquad) * conv_num_flux_dot_n[iquad][istate] * JxW[iquad];
+                    // Diffusive
+                    rhs[istate] -= fe_values_boundary->shape_value(itest,iquad) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW[iquad];
+                    rhs[istate] += fe_values_boundary->shape_grad(itest,iquad) * diss_flux_jump_int[iquad][istate] * JxW[iquad];
                 }
+                //for (int istate=0; istate<nstate; istate++) {
+                //    rhs[istate] -= 
+                //        diss_soln_num_flux[iquad][istate] *
+                //        fe_values_boundary->shape_value(itest,iquad) *
+                //        JxW[iquad];
+                //}
 
                 //     // Diffusion
                 //     // *******************
@@ -219,7 +256,7 @@ namespace PHiLiP
                 //     // Weakly imposed boundary condition
                 //     //soln_ext = boundary_values[iquad];
 
-                //     const ADArrayVector soln_grad_ext = soln_grad_int[iquad];
+                //     const ADArrayTensor1 soln_grad_ext = soln_grad_int[iquad];
 
                 //     // Note: The test function is piece-wise defined to be non-zero only on the associated cell
                 //     //       and zero everywhere else.
@@ -229,9 +266,9 @@ namespace PHiLiP
                 //     const Tensor<1,dim,real> test_grad_ext;// = 0;
 
                 //     //const Tensor<1,dim,ADtype> soln_jump        = (soln_int[iquad] - soln_ext) * normal_int;
-                //     const ADArrayVector soln_grad_avg    = 0.5*(soln_grad_int[iquad] + soln_grad_ext);
-                //     const ADArrayVector test_jump        = (test_int - test_ext) * normal_int;
-                //     const ADArrayVector test_grad_avg    = 0.5*(test_grad_int + test_grad_ext);
+                //     const ADArrayTensor1 soln_grad_avg    = 0.5*(soln_grad_int[iquad] + soln_grad_ext);
+                //     const ADArrayTensor1 test_jump        = (test_int - test_ext) * normal_int;
+                //     const ADArrayTensor1 test_grad_avg    = 0.5*(test_grad_int + test_grad_ext);
 
                 //     //rhs += (soln_jump * test_grad_avg + soln_grad_avg * test_jump - penalty*soln_jump*test_jump) * JxW[iquad];
             }
@@ -265,7 +302,7 @@ namespace PHiLiP
     {
         using ADtype = Sacado::Fad::DFad<real>;
         using ADArray = std::array<ADtype,nstate>;
-        using ADArrayVector = std::array< Tensor<1,dim,ADtype>, nstate >;
+        using ADArrayTensor1 = std::array< Tensor<1,dim,ADtype>, nstate >;
 
         // Use quadrature points of neighbor cell
         // Might want to use the maximum n_quad_pts1 and n_quad_pts2
@@ -304,18 +341,26 @@ namespace PHiLiP
         std::vector<real> dR2_dW1(n_dofs_current_cell);
         std::vector<real> dR2_dW2(n_dofs_neighbor_cell);
 
-        std::vector<ADArray> normal_int_numerical_flux(n_face_quad_pts);
+        std::vector<ADArray> conv_num_flux_dot_n(n_face_quad_pts);
 
         // Interpolate solution to the face quadrature points
         std::vector< ADArray > soln_int(n_face_quad_pts);
         std::vector< ADArray > soln_ext(n_face_quad_pts);
 
-        std::vector< ADArrayVector > soln_grad_int(n_face_quad_pts); // Tensor initialize with zeros
-        std::vector< ADArrayVector > soln_grad_ext(n_face_quad_pts); // Tensor initialize with zeros
+        std::vector< ADArrayTensor1 > soln_grad_int(n_face_quad_pts); // Tensor initialize with zeros
+        std::vector< ADArrayTensor1 > soln_grad_ext(n_face_quad_pts); // Tensor initialize with zeros
+
+        std::vector<ADArray> diss_soln_num_flux(n_face_quad_pts); // u*
+        std::vector<ADArray> diss_auxi_num_flux_dot_n(n_face_quad_pts); // sigma*
+
+        std::vector<ADArrayTensor1> diss_flux_jump_int(n_face_quad_pts); // u*-u_int
+        std::vector<ADArrayTensor1> diss_flux_jump_ext(n_face_quad_pts); // u*-u_ext
 
         for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
 
             const Tensor<1,dim,ADtype> normal_int = normals_int[iquad];
+            const Tensor<1,dim,ADtype> normal_ext = -normal_int;
+
 
             for (int istate=0; istate<nstate; ++istate) {
                 soln_int[iquad][istate]      = 0;
@@ -333,11 +378,26 @@ namespace PHiLiP
                     soln_grad_ext[iquad][istate] += neighbor_solution_ad[itrial][istate] * fe_values_ext->shape_grad(itrial, iquad);
                 }
             }
-            normal_int_numerical_flux[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+            conv_num_flux_dot_n[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+
+            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;
+            for (int s=0; s<nstate; s++) {
+                diss_soln_jump_int[s] = (diss_soln_num_flux[iquad][s] - soln_int[iquad][s]) * normal_int;
+                diss_soln_jump_ext[s] = (diss_soln_num_flux[iquad][s] - soln_ext[iquad][s]) * normal_ext;
+            }
+            pde_physics->dissipative_flux (soln_int[iquad], diss_soln_jump_int, diss_flux_jump_int[iquad]);
+            pde_physics->dissipative_flux (soln_ext[iquad], diss_soln_jump_ext, diss_flux_jump_ext[iquad]);
+
+            diss_auxi_num_flux_dot_n[iquad] = diss_num_flux->evaluate_auxiliary_flux(
+                soln_int[iquad], soln_ext[iquad],
+                soln_grad_int[iquad], soln_grad_ext[iquad],
+                normal_int, penalty);
 
         }
 
-        for (unsigned int itest_current=0; itest_current<n_dofs_current_cell; ++itest_current) {
+        for (unsigned int itest_int=0; itest_int<n_dofs_current_cell; ++itest_int) {
             // From test functions associated with current cell point of view
             // *******************
             ADArray rhs;
@@ -346,67 +406,29 @@ namespace PHiLiP
             }
 
             for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
-                // Convection
                 for (int istate=0; istate<nstate; istate++) {
-                    rhs[istate] -=
-                        fe_values_int->shape_value(itest_current,iquad) *
-                        normal_int_numerical_flux[iquad][istate] *
-                        JxW_int[iquad];
+                    // Convection
+                    rhs[istate] -= fe_values_int->shape_value(itest_int,iquad) * conv_num_flux_dot_n[iquad][istate] * JxW_int[iquad];
+                    // Diffusive
+                    rhs[istate] -= fe_values_int->shape_value(itest_int,iquad) * diss_auxi_num_flux_dot_n[iquad][istate] * JxW_int[iquad];
+                    rhs[istate] += fe_values_int->shape_grad(itest_int,iquad) * diss_flux_jump_int[iquad][istate] * JxW_int[iquad];
                 }
-
-                //    // Diffusion
-
-
-                //    //// Note: The test function is piece-wise defined to be non-zero only on the associated cell
-                //    ////       and zero everywhere else.
-                //    //const Tensor< 1, dim, real > test_grad_int = fe_values_int->shape_grad(itest_current, iquad);
-                //    //const Tensor< 1, dim, real > test_grad_ext; // Constructor initializes with zeros
-
-                //    //const Tensor< 1, dim, real > normal1        = normals_int[iquad];
-                //    //const Tensor< 1, dim, ADtype > soln_jump      = (soln_int[iquad] - soln_ext[iquad]) * normal1;
-                //    //const Tensor< 1, dim, ADtype > soln_grad_avg  = 0.5*(soln_grad_ext[iquad] + soln_grad_int[iquad]);
-
-                //    //const Tensor< 1, dim, real > test_jump        = (test_int - test_ext) * normal1;
-                //    //const Tensor< 1, dim, real > test_grad_avg    = 0.5*(test_grad_int + test_grad_ext);
-
-                //    //rhs += (soln_jump * test_grad_avg + soln_grad_avg * test_jump - penalty*soln_jump*test_jump) * JxW_int[iquad];
-
-                //    // Note: The test function is piece-wise defined to be non-zero only on the associated cell
-                //    //       and zero everywhere else.
-                //    const real test_int = fe_values_int->shape_value(itest_current, iquad);
-                //    const real test_ext = 0;
-
-                //    for (unsigned int idim=0; idim<dim; ++idim) {
-                //        // Note: The test function is piece-wise defined to be non-zero only on the associated cell
-                //        //       and zero everywhere else.
-                //        const real test_grad_int = fe_values_int->shape_grad(itest_current, iquad)[idim];
-                //        const real test_grad_ext = 0;
-
-                //        // Using normals from interior point of view
-                //        const real normal1          = normals_int[iquad][idim];
-                //        const ADtype soln_jump      = (soln_int[iquad] - soln_ext[iquad]) * normal1;
-                //        const ADtype soln_grad_avg  = 0.5*(soln_grad_int[iquad][idim] + soln_grad_ext[iquad][idim]);
-                //        const real test_jump        = (test_int - test_ext) * normal1;
-                //        const real test_grad_avg    = 0.5*(test_grad_int + test_grad_ext);
-                //    
-                //        //rhs += (soln_jump * test_grad_avg + soln_grad_avg * test_jump - penalty*soln_jump*test_jump) * JxW_int[iquad];
-                //    }
             }
             // *******************
 
             const int ISTATE = 0;
-            current_cell_rhs(itest_current) += rhs[ISTATE].val();
+            current_cell_rhs(itest_int) += rhs[ISTATE].val();
             for (unsigned int itrial = 0; itrial < n_dofs_current_cell; ++itrial) {
                 dR1_dW1[itrial] = rhs[ISTATE].dx(itrial);
             }
             for (unsigned int itrial = 0; itrial < n_dofs_neighbor_cell; ++itrial) {
                 dR1_dW2[itrial] = rhs[ISTATE].dx(n_dofs_current_cell+itrial);
             }
-            system_matrix.add(current_dofs_indices[itest_current], current_dofs_indices, dR1_dW1);
-            system_matrix.add(current_dofs_indices[itest_current], neighbor_dofs_indices, dR1_dW2);
+            system_matrix.add(current_dofs_indices[itest_int], current_dofs_indices, dR1_dW1);
+            system_matrix.add(current_dofs_indices[itest_int], neighbor_dofs_indices, dR1_dW2);
         }
 
-        for (unsigned int itest_neighbor=0; itest_neighbor<n_dofs_neighbor_cell; ++itest_neighbor) {
+        for (unsigned int itest_ext=0; itest_ext<n_dofs_neighbor_cell; ++itest_ext) {
             // From test functions associated with neighbour cell point of view
             // *******************
             ADArray rhs;
@@ -415,50 +437,25 @@ namespace PHiLiP
             }
 
             for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
-                // Convection
                 for (int istate=0; istate<nstate; istate++) {
-                    rhs[istate] -=
-                        fe_values_ext->shape_value(itest_neighbor,iquad) *
-                        (-normal_int_numerical_flux[iquad][istate]) *
-                        JxW_int[iquad];
+                    // Convection
+                    rhs[istate] -= fe_values_ext->shape_value(itest_ext,iquad) * (-conv_num_flux_dot_n[iquad][istate]) * JxW_int[iquad];
+                    // Diffusive
+                    rhs[istate] -= fe_values_ext->shape_value(itest_ext,iquad) * (-diss_auxi_num_flux_dot_n[iquad][istate]) * JxW_int[iquad];
+                    rhs[istate] += fe_values_ext->shape_grad(itest_ext,iquad) * diss_flux_jump_ext[iquad][istate] * JxW_int[iquad];
                 }
-
-                //    // Diffusion
-
-                //    // Note: The test function is piece-wise defined to be non-zero only on the associated cell
-                //    //       and zero everywhere else.
-                //    const real test_int = fe_values_ext->shape_value(itest_neighbor, iquad);
-                //    const real test_ext = 0;
-
-                //    for (unsigned int idim=0; idim<dim; ++idim) {
-                //        // Note: The test function is piece-wise defined to be non-zero only on the associated cell
-                //        //       and zero everywhere else.
-                //        const real test_grad_int = fe_values_ext->shape_grad(itest_neighbor, iquad)[idim];
-                //        const real test_grad_ext = 0;
-
-                //        // Using normals from other side
-                //        // So opposite and equal value
-                //        const real flipped_normal1 = -normals_int[iquad][idim];
-                //        // Flipped soln_ext and soln_int
-                //        const ADtype soln_jump      = (soln_ext[iquad] - soln_int[iquad]) * flipped_normal1;
-                //        const ADtype soln_grad_avg  = 0.5*(soln_grad_int[iquad][idim] + soln_grad_ext[iquad][idim]);
-                //        const real test_jump        = (test_int - test_ext) * flipped_normal1;
-                //        const real test_grad_avg    = 0.5*(test_grad_int + test_grad_ext);
-                //    
-                //        //rhs += (soln_jump * test_grad_avg + soln_grad_avg * test_jump - penalty*soln_jump*test_jump) * JxW_int[iquad];
-                //    }
             }
             // *******************
             const int ISTATE = 0;
-            neighbor_cell_rhs(itest_neighbor) += rhs[ISTATE].val();
+            neighbor_cell_rhs(itest_ext) += rhs[ISTATE].val();
             for (unsigned int itrial = 0; itrial < n_dofs_current_cell; ++itrial) {
                 dR2_dW1[itrial] = rhs[ISTATE].dx(itrial);
             }
             for (unsigned int itrial = 0; itrial < n_dofs_neighbor_cell; ++itrial) {
                 dR2_dW2[itrial] = rhs[ISTATE].dx(n_dofs_current_cell+itrial);
             }
-            system_matrix.add(neighbor_dofs_indices[itest_neighbor], current_dofs_indices, dR2_dW1);
-            system_matrix.add(neighbor_dofs_indices[itest_neighbor], neighbor_dofs_indices, dR2_dW2);
+            system_matrix.add(neighbor_dofs_indices[itest_ext], current_dofs_indices, dR2_dW1);
+            system_matrix.add(neighbor_dofs_indices[itest_ext], neighbor_dofs_indices, dR2_dW2);
         }
     }
     template void DiscontinuousGalerkin<PHILIP_DIM, 1, double>::assemble_face_term_implicit(
diff --git a/src/dg.cpp b/src/dg.cpp
index 503dfd85c..66752656a 100644
--- a/src/dg.cpp
+++ b/src/dg.cpp
@@ -70,6 +70,8 @@ namespace PHiLiP
         pde_physics = PhysicsFactory<dim, nstate, ADtype >::create_Physics(parameters->pde_type);
         conv_num_flux = NumericalFluxFactory<dim, nstate, ADtype>::create_convective_numerical_flux
             (parameters->conv_num_flux_type, pde_physics);
+        diss_num_flux = NumericalFluxFactory<dim, nstate, ADtype>::create_dissipative_numerical_flux
+            (parameters->diss_num_flux_type, pde_physics);
     }
     // Destructor
     template <int dim, int nstate, typename real>
@@ -99,6 +101,7 @@ namespace PHiLiP
     void DiscontinuousGalerkin<dim,nstate,real>::delete_fe_values ()
     {
         delete conv_num_flux;
+        delete diss_num_flux;
         delete pde_physics;
 
         if (fe_values_cell          != NULL) delete fe_values_cell;
@@ -538,6 +541,44 @@ namespace PHiLiP
       data_out.write_gnuplot(gnuplot_output);
     }
 
+    template <int dim, int nstate, typename real>
+    void DiscontinuousGalerkin<dim,nstate,real>::evaluate_inverse_mass_matrices ()
+    {
+        // Invert and store mass matrix
+        // Using Gauss-Jordan since it's deal.II's default invert function
+        // Could store Cholesky decomposition for more efficient pre-processing
+        const int n_quad_pts      = quadrature.size();
+        const int n_dofs_per_cell = fe.dofs_per_cell;
+
+        typename DoFHandler<dim>::active_cell_iterator
+           cell = dof_handler.begin_active(),
+           endc = dof_handler.end();
+
+        inv_mass_matrix.resize(triangulation->n_active_cells(),
+                               FullMatrix<real>(n_dofs_per_cell));
+        FullMatrix<real> mass_matrix(n_dofs_per_cell);
+        for (; cell!=endc; ++cell) {
+
+            int cell_index = cell->index();
+            fe_values_cell->reinit(cell);
+
+            for (int itest=0; itest<n_dofs_per_cell; ++itest) {
+                for (int itrial=itest; itrial<n_dofs_per_cell; ++itrial) {
+                    mass_matrix[itest][itrial] = 0.0;
+                    for (int iquad=0; iquad<n_quad_pts; ++iquad) {
+                        mass_matrix[itest][itrial] +=
+                            fe_values_cell->shape_value(itest,iquad)
+                            * fe_values_cell->shape_value(itrial,iquad)
+                            * fe_values_cell->JxW(iquad);
+                    }
+                    mass_matrix[itrial][itest] = mass_matrix[itest][itrial];
+                }
+            }
+            inv_mass_matrix[cell_index].invert(mass_matrix);
+        }
+    }
+
+
 
 
     template class DiscontinuousGalerkin <PHILIP_DIM, 1, double>;
diff --git a/src/dg.h b/src/dg.h
index c056508a9..4b4411dc5 100644
--- a/src/dg.h
+++ b/src/dg.h
@@ -44,6 +44,7 @@ namespace PHiLiP
 
         void allocate_system ();
         void assemble_system ();
+        void evaluate_inverse_mass_matrices ();
         double get_residual_l2norm ();
 
         void delete_fe_values ();
@@ -84,6 +85,7 @@ namespace PHiLiP
         /// Contains the physics of the PDE
         Physics<dim, nstate, Sacado::Fad::DFad<real> >  *pde_physics;
         NumericalFluxConvective<dim, nstate, Sacado::Fad::DFad<real> > *conv_num_flux;
+        NumericalFluxDissipative<dim, nstate, Sacado::Fad::DFad<real> > *diss_num_flux;
 
 
         void allocate_system_explicit ();
@@ -141,6 +143,8 @@ namespace PHiLiP
             Vector<real>          &current_cell_rhs,
             Vector<real>          &neighbor_cell_rhs);
 
+        std::vector<FullMatrix<real>> inv_mass_matrix;
+
         // QGauss is Gauss-Legendre quadrature nodes
         const QGauss<dim>   quadrature;
         const QGauss<dim-1> face_quadrature;
diff --git a/src/grid_study.cpp b/src/grid_study.cpp
index bfcd819f5..35f69319e 100644
--- a/src/grid_study.cpp
+++ b/src/grid_study.cpp
@@ -90,7 +90,7 @@ namespace PHiLiP
 
             // p0 tends to require a finer grid to reach asymptotic region
             unsigned int n_grids = n_grids_input;
-            if (poly_degree == 0) n_grids = n_grids_input + 2;
+            if (poly_degree <= 1) n_grids = n_grids_input + 1;
 
             std::vector<int> n_1d_cells(n_grids);
             n_1d_cells[0] = initial_grid_size;
@@ -168,6 +168,8 @@ namespace PHiLiP
                 //dg.set_physics(physics);
                 dg.allocate_system ();
 
+                dg.evaluate_inverse_mass_matrices();
+
                 //ODESolver<dim, double> *ode_solver = ODESolverFactory<dim, double>::create_ODESolver(parameters.solver_type);
                 ODESolver<dim, nstate, double> *ode_solver = ODESolverFactory<dim, nstate, double>::create_ODESolver(&dg);
 
diff --git a/src/linear_solver.cpp b/src/linear_solver.cpp
index fc45a14cc..c61d96c09 100644
--- a/src/linear_solver.cpp
+++ b/src/linear_solver.cpp
@@ -20,15 +20,18 @@ namespace PHiLiP
 
         //system_matrix.print(std::cout, true);
         //solution.print(std::cout);
-        //right_hand_side.print(std::cout);
-        //FullMatrix<double> fullA(system_matrix.m());
-        //fullA.copy_from(system_matrix);
-        //std::cout<<"Dense matrix:"<<std::endl;
-        //fullA.print_formatted(std::cout, 3, true, 10, "0", 1., 0.);
 
 
         using LinSolvParam = Parameters::LinearSolver;
         if (param->linear_solver_type == LinSolvParam::LinearSolverType::direct) {
+            if (param->linear_solver_output == Parameters::OutputType::verbose) {
+                right_hand_side.print(std::cout);
+                FullMatrix<double> fullA(system_matrix.m());
+                fullA.copy_from(system_matrix);
+                std::cout<<"Dense matrix:"<<std::endl;
+                fullA.print_formatted(std::cout, 3, true, 10, "0", 1., 0.);
+            }
+
             SolverControl solver_control(1, 0);
             TrilinosWrappers::SolverDirect::AdditionalData data(false);
             //TrilinosWrappers::SolverDirect::AdditionalData data(parameters.output == Parameters::Solver::verbose);
diff --git a/src/numerical_flux/CMakeLists.txt b/src/numerical_flux/CMakeLists.txt
index 8e77b0498..1a089b422 100644
--- a/src/numerical_flux/CMakeLists.txt
+++ b/src/numerical_flux/CMakeLists.txt
@@ -1,4 +1,5 @@
 set(SOURCE
+    viscous_numerical_flux.cpp
     numerical_flux.cpp
     )
 
diff --git a/src/numerical_flux/numerical_flux.cpp b/src/numerical_flux/numerical_flux.cpp
index 20de4f0ac..6bebdae8f 100644
--- a/src/numerical_flux/numerical_flux.cpp
+++ b/src/numerical_flux/numerical_flux.cpp
@@ -1,6 +1,7 @@
 #include <Sacado.hpp>
 #include <deal.II/differentiation/ad/sacado_product_types.h>
 #include "numerical_flux.h"
+#include "viscous_numerical_flux.h"
 
 namespace PHiLiP
 {
@@ -87,30 +88,6 @@ namespace PHiLiP
         return numerical_flux_dot_n;
     }
 
-    template <int dim, int nstate, typename real>
-    NumericalFluxDissipative<dim,nstate,real>::~NumericalFluxDissipative() {}
-
-    template<int dim, int nstate, typename real>
-    void SymmetricInternalPenalty<dim,nstate,real>
-    ::evaluate_flux (
-        const std::array<real, nstate> &soln_int,
-        const std::array<real, nstate> &soln_ext,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
-        const Tensor<1,dim,real> &normal1,
-        std::array<real, nstate> &soln_flux,
-        std::array<Tensor<1,dim,real>, nstate> &grad_flux) const
-    {
-        std::array<real,nstate> soln_avg = array_average<nstate,real>(soln_int, soln_ext);
-
-        //std::array<Tensor<2,dim,real>,nstate> diffusion_matrix_int = 
-        //    pde_physics->diffusion_matrix(soln_int);
-
-        //std::array<Tensor<2,dim,real>,nstate> diffusion_matrix_ext = 
-        //    pde_physics->diffusion_matrix(soln_ext);
-    }
-
-
     // Instantiation
     template class NumericalFluxConvective<PHILIP_DIM, 1, double>;
     template class NumericalFluxConvective<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
@@ -118,11 +95,6 @@ namespace PHiLiP
     template class LaxFriedrichs<PHILIP_DIM, 1, double>;
     template class LaxFriedrichs<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
 
-    template class NumericalFluxDissipative<PHILIP_DIM, 1, double>;
-    template class NumericalFluxDissipative<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
-
-    template class SymmetricInternalPenalty<PHILIP_DIM, 1, double>;
-    template class SymmetricInternalPenalty<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
 
     template class NumericalFluxFactory<PHILIP_DIM, 1, double>;
     template class NumericalFluxFactory<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
diff --git a/src/numerical_flux/numerical_flux.h b/src/numerical_flux/numerical_flux.h
index c797c59fb..08d85456b 100644
--- a/src/numerical_flux/numerical_flux.h
+++ b/src/numerical_flux/numerical_flux.h
@@ -3,6 +3,7 @@
 
 #include <deal.II/base/tensor.h>
 #include "physics/physics.h"
+#include "numerical_flux/viscous_numerical_flux.h"
 namespace PHiLiP
 {
     using namespace dealii;
@@ -46,56 +47,6 @@ namespace PHiLiP
 
     };
 
-    /// Numerical flux associated with dissipation
-    template<int dim, int nstate, typename real>
-    class NumericalFluxDissipative
-    {
-    public:
-    /// Base class destructor.
-    /// Abstract classes required virtual destructors.
-    /// Even if it is a pure function, its definition is still required.
-    virtual ~NumericalFluxDissipative() = 0;
-
-    virtual void evaluate_flux (
-        const std::array<real, nstate> &soln_int,
-        const std::array<real, nstate> &soln_ext,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
-        const Tensor<1,dim,real> &normal1,
-        std::array<real, nstate> &soln_flux,
-        std::array<Tensor<1,dim,real>, nstate> &grad_flux) const = 0;
-
-    };
-
-    template<int dim, int nstate, typename real>
-    class SymmetricInternalPenalty: public NumericalFluxDissipative<dim, nstate, real>
-    {
-    public:
-    /// Constructor
-    SymmetricInternalPenalty(Physics<dim, nstate, real> *physics_input)
-    :
-    pde_physics(physics_input)
-    {};
-    /// Destructor
-    ~SymmetricInternalPenalty() {};
-
-    /// Evaluate solution and gradient flux
-    /*  $\hat{u} = {u_h}$, 
-     *  $ \hat{A} = {{ A \nabla u_h }} - \mu {{ A }} [[ u_h ]] $
-     */
-    void evaluate_flux (
-        const std::array<real, nstate> &soln_int,
-        const std::array<real, nstate> &soln_ext,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
-        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
-        const Tensor<1,dim,real> &normal1,
-        std::array<real, nstate> &soln_flux,
-        std::array<Tensor<1,dim,real>, nstate> &grad_flux) const;
-    protected:
-    const Physics<dim, nstate, real> *pde_physics;
-
-    };
-
 
     template <int dim, int nstate, typename real>
     class NumericalFluxFactory
diff --git a/src/numerical_flux/viscous_numerical_flux.cpp b/src/numerical_flux/viscous_numerical_flux.cpp
new file mode 100644
index 000000000..7f1ac5b64
--- /dev/null
+++ b/src/numerical_flux/viscous_numerical_flux.cpp
@@ -0,0 +1,157 @@
+#include <Sacado.hpp>
+#include <deal.II/differentiation/ad/sacado_product_types.h>
+#include "viscous_numerical_flux.h"
+
+namespace PHiLiP
+{
+    using namespace dealii;
+    using AllParam = Parameters::AllParameters;
+
+    // Protyping low level functions
+    template<int nstate, typename real_tensor>
+    std::array<real_tensor, nstate> array_average(
+        const std::array<real_tensor, nstate> &array1,
+        const std::array<real_tensor, nstate> &array2)
+    {
+        std::array<real_tensor,nstate> array_average;
+        for (int s=0; s<nstate; s++) {
+            array_average[s] = 0.5*(array1[s] + array2[s]);
+        }
+        return array_average;
+    }
+
+    template<int dim, int nstate, typename real>
+    std::array<Tensor<1,dim,real>,nstate> array_jump(
+        const std::array<real, nstate> &array1,
+        const std::array<real, nstate> &array2,
+        const Tensor<1,dim,real> &normal1)
+    {
+        std::array<Tensor<1,dim,real>,nstate> array_jump;
+        for (int s=0; s<nstate; s++) {
+            array_jump[s] = (array1[s] - array2[s])*normal1;
+        }
+        return array_jump;
+    }
+
+    template <int dim, int nstate, typename real>
+    NumericalFluxDissipative<dim,nstate,real>::~NumericalFluxDissipative() {}
+
+
+    //template<int dim, int nstate, typename real>
+    //void NumericalFluxDissipative<dim,nstate,real>
+    //::evaluate_auxiliary (
+    //    const std::vector<std::array<real, nstate>> &soln_nodal_face_int,
+    //    const std::vector<std::array<real, nstate>> &soln_nodal_face_ext,
+    //    const std::vector<std::array<real, nstate>> &soln_nodal_face_flux,
+    //    const std::vector<std::array<real, nstate>> &soln_grad_nodal_volume_int,
+    //    const std::vector<std::array<real, nstate>> &soln_grad_nodal_volume_ext,
+    //    const std::vector<Tensor<1,dim,real>> &normal_int,
+    //    std::vector<std::array<real, nstate>> &auxiliary_modal_int,
+    //    std::vector<std::array<real, nstate>> &auxiliary_modal_ext) const
+    //{
+    //    // Need some kind of assert to confirm that A and A_transpose correspond
+    //    // to the current interior and exterior cell
+
+    //    // Interpolate solution to the face quadrature points
+    //    std::vector< ADArray > soln_int(n_face_quad_pts);
+    //    std::vector< ADArray > soln_ext(n_face_quad_pts);
+
+    //    std::vector< ADArrayVector > soln_grad_int(n_face_quad_pts); // Tensor initialize with zeros
+    //    std::vector< ADArrayVector > soln_grad_ext(n_face_quad_pts); // Tensor initialize with zeros
+    //    for (unsigned int iquad=0; iquad<n_face_quad_pts; ++iquad) {
+
+    //        const Tensor<1,dim,ADtype> normal_int = normals_int[iquad];
+
+    //        for (int istate=0; istate<nstate; ++istate) {
+    //            soln_int[iquad][istate]      = 0;
+    //            soln_ext[iquad][istate]      = 0;
+    //            soln_grad_int[iquad][istate] = 0;
+    //            soln_grad_ext[iquad][istate] = 0;
+
+    //            // Interpolate solution to face
+    //            for (unsigned int itrial=0; itrial<n_dofs_current_cell; itrial++) {
+    //                soln_int[iquad][istate]      += current_solution_ad[itrial][istate] * fe_values_int->shape_value(itrial, iquad);
+    //                soln_grad_int[iquad][istate] += current_solution_ad[itrial][istate] * fe_values_int->shape_grad(itrial, iquad);
+    //            }
+    //            for (unsigned int itrial=0; itrial<n_dofs_neighbor_cell; itrial++) {
+    //                soln_ext[iquad][istate]      += neighbor_solution_ad[itrial][istate] * fe_values_ext->shape_value(itrial, iquad);
+    //                soln_grad_ext[iquad][istate] += neighbor_solution_ad[itrial][istate] * fe_values_ext->shape_grad(itrial, iquad);
+    //            }
+    //        }
+    //        normal_int_numerical_flux[iquad] = conv_num_flux->evaluate_flux(soln_int[iquad], soln_ext[iquad], normal_int);
+
+    //    }
+    //}
+
+
+    //template<int dim, int nstate, typename real>
+    //void NumericalFluxDissipative<dim,nstate,real>
+    //::assemble_auxiliary_cell (
+    //    const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int) const
+    //{
+    //}
+
+    template<int dim, int nstate, typename real>
+    std::array<real, nstate> SymmetricInternalPenalty<dim,nstate,real>
+    ::evaluate_solution_flux (
+        const std::array<real, nstate> &soln_int,
+        const std::array<real, nstate> &soln_ext,
+        const Tensor<1,dim,real> &normal_int) const
+    {
+        std::array<real,nstate> soln_avg = array_average<nstate,real>(soln_int, soln_ext);
+
+        //std::array<Tensor<2,dim,real>,nstate> diffusion_matrix_int = 
+        //    pde_physics->diffusion_matrix(soln_int);
+
+        //std::array<Tensor<2,dim,real>,nstate> diffusion_matrix_ext = 
+        //    pde_physics->diffusion_matrix(soln_ext);
+
+        return soln_avg;
+    }
+
+    template<int dim, int nstate, typename real>
+    std::array<real, nstate> SymmetricInternalPenalty<dim,nstate,real>
+    ::evaluate_auxiliary_flux (
+        const std::array<real, nstate> &soln_int,
+        const std::array<real, nstate> &soln_ext,
+        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
+        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
+        const Tensor<1,dim,real> &normal_int,
+        const real &penalty) const
+    {
+        using ArrayTensor1 = std::array<Tensor<1,dim,real>, nstate>;
+        ArrayTensor1 phys_flux_int, phys_flux_ext;
+
+        // {{A*grad_u}}
+        pde_physics->dissipative_flux (soln_int, soln_grad_int, phys_flux_int);
+        pde_physics->dissipative_flux (soln_ext, soln_grad_ext, phys_flux_ext);
+        ArrayTensor1 phys_flux_avg = array_average<nstate,Tensor<1,dim,real>>(phys_flux_int, phys_flux_ext);
+
+        // {{A}}*[[u]]
+        ArrayTensor1 soln_jump     = array_jump<dim,nstate,real>(soln_int, soln_ext, normal_int);
+        ArrayTensor1 A_jumpu_int, A_jumpu_ext;
+        pde_physics->dissipative_flux (soln_int, soln_jump, A_jumpu_int);
+        pde_physics->dissipative_flux (soln_ext, soln_jump, A_jumpu_ext);
+        const ArrayTensor1 A_jumpu_avg = array_average<nstate,Tensor<1,dim,real>>(A_jumpu_int, A_jumpu_ext);
+
+
+        std::array<real,nstate> auxiliary_flux_dot_n;
+        for (int s=0; s<nstate; s++) {
+            auxiliary_flux_dot_n[s] = (phys_flux_avg[s] - penalty * A_jumpu_avg[s]) * normal_int;
+            //auxiliary_flux_dot_n[s] = (phys_flux_avg[s] - penalty * soln_jump[s]) * normal_int;
+        }
+        return auxiliary_flux_dot_n;
+
+    }
+
+
+
+
+    // Instantiation
+    template class NumericalFluxDissipative<PHILIP_DIM, 1, double>;
+    template class NumericalFluxDissipative<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
+
+    template class SymmetricInternalPenalty<PHILIP_DIM, 1, double>;
+    template class SymmetricInternalPenalty<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
+
+}
diff --git a/src/numerical_flux/viscous_numerical_flux.h b/src/numerical_flux/viscous_numerical_flux.h
new file mode 100644
index 000000000..8e27dd606
--- /dev/null
+++ b/src/numerical_flux/viscous_numerical_flux.h
@@ -0,0 +1,77 @@
+#ifndef __VISCOUS_NUMERICAL_FLUX__
+#define __VISCOUS_NUMERICAL_FLUX__
+
+#include <deal.II/base/tensor.h>
+#include "physics/physics.h"
+namespace PHiLiP
+{
+    using namespace dealii;
+    using AllParam = Parameters::AllParameters;
+
+    /// Numerical flux associated with dissipation
+    template<int dim, int nstate, typename real>
+    class NumericalFluxDissipative
+    {
+    public:
+    /// Base class destructor.
+    /// Abstract classes required virtual destructors.
+    /// Even if it is a pure function, its definition is still required.
+    virtual ~NumericalFluxDissipative() = 0;
+
+    virtual std::array<real, nstate> evaluate_solution_flux (
+        const std::array<real, nstate> &soln_int,
+        const std::array<real, nstate> &soln_ext,
+        const Tensor<1,dim,real> &normal_int) const = 0;
+
+    virtual std::array<real, nstate> evaluate_auxiliary_flux (
+        const std::array<real, nstate> &soln_int,
+        const std::array<real, nstate> &soln_ext,
+        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
+        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
+        const Tensor<1,dim,real> &normal_int,
+        const real &penalty) const = 0;
+
+    Tensor<1,dim, Tensor<1,nstate,real>> diffusion_matrix_int;
+    Tensor<1,dim, Tensor<1,nstate,real>> diffusion_matrix_int_transpose;
+    Tensor<1,dim, Tensor<1,nstate,real>> diffusion_matrix_ext;
+    Tensor<1,dim, Tensor<1,nstate,real>> diffusion_matrix_ext_transpose;
+
+    };
+
+    template<int dim, int nstate, typename real>
+    class SymmetricInternalPenalty: public NumericalFluxDissipative<dim, nstate, real>
+    {
+    public:
+    /// Constructor
+    SymmetricInternalPenalty(Physics<dim, nstate, real> *physics_input)
+    :
+    pde_physics(physics_input)
+    {};
+    /// Destructor
+    ~SymmetricInternalPenalty() {};
+
+    /// Evaluate solution and gradient flux
+    /*  $\hat{u} = {u_h}$, 
+     *  $ \hat{A} = {{ A \nabla u_h }} - \mu {{ A }} [[ u_h ]] $
+     */
+    std::array<real, nstate> evaluate_solution_flux (
+        const std::array<real, nstate> &soln_int,
+        const std::array<real, nstate> &soln_ext,
+        const Tensor<1,dim,real> &normal_int) const;
+
+    std::array<real, nstate> evaluate_auxiliary_flux (
+        const std::array<real, nstate> &soln_int,
+        const std::array<real, nstate> &soln_ext,
+        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_int,
+        const std::array<Tensor<1,dim,real>, nstate> &soln_grad_ext,
+        const Tensor<1,dim,real> &normal_int,
+        const real &penalty) const;
+        
+    protected:
+    const Physics<dim, nstate, real> *pde_physics;
+
+    };
+
+}
+
+#endif
diff --git a/src/ode_solver.cpp b/src/ode_solver.cpp
index 62874883b..3a8ef7caa 100644
--- a/src/ode_solver.cpp
+++ b/src/ode_solver.cpp
@@ -16,25 +16,39 @@ namespace PHiLiP
         dg->assemble_system ();
         this->residual_norm = dg->get_residual_l2norm();
         this->residual_norm = 1; // Always do at least 1 iteration
+        double update_norm = 1; // Always do at least 1 iteration
         this->current_iteration = 0;
 
-        while (    this->residual_norm     > parameters->nonlinear_steady_residual_tolerance 
+        while (    
+                   this->residual_norm     > parameters->nonlinear_steady_residual_tolerance 
+                && update_norm             > parameters->nonlinear_steady_residual_tolerance 
                 && this->current_iteration < parameters->nonlinear_max_iterations )
         {
-            ++this->current_iteration;
-
             if ((this->parameters->ode_output) == Parameters::OutputType::verbose &&
                 (this->current_iteration%parameters->print_iteration_modulo) == 0 )
             std::cout << " Iteration: " << this->current_iteration 
                       << " Residual norm: " << this->residual_norm
                       << std::endl;
 
+            if ((this->parameters->ode_output) == Parameters::OutputType::verbose &&
+                (this->current_iteration%parameters->print_iteration_modulo) == 0 )
+            std::cout << " Assembling system... ";
+            dg->assemble_system ();
+
+            if ((this->parameters->ode_output) == Parameters::OutputType::verbose &&
+                (this->current_iteration%parameters->print_iteration_modulo) == 0 )
+            std::cout << " Evaluating system update... ";
             evaluate_solution_update ();
+
+
             //this->solution_update *= 0.1;
+            update_norm = this->solution_update.l2_norm();
             dg->solution += this->solution_update;
 
-            dg->assemble_system ();
             this->residual_norm = dg->get_residual_l2norm();
+
+            ++(this->current_iteration);
+
         }
         return 1;
     }
diff --git a/src/parameters.cpp b/src/parameters.cpp
index 255f4e38a..f1106fb54 100644
--- a/src/parameters.cpp
+++ b/src/parameters.cpp
@@ -74,7 +74,7 @@ namespace Parameters
     {
         prm.enter_subsection("ODE solver");
         {
-            prm.declare_entry("ode_output", "quiet",
+            prm.declare_entry("ode_output", "verbose",
                               Patterns::Selection("quiet|verbose"),
                               "State whether output from ODE solver should be printed. "
                               "Choices are <quiet|verbose>.");
diff --git a/src/physics/CMakeLists.txt b/src/physics/CMakeLists.txt
index ef3875755..b828b426c 100644
--- a/src/physics/CMakeLists.txt
+++ b/src/physics/CMakeLists.txt
@@ -1,7 +1,7 @@
 set(SOURCE
-    linear_advection.cpp
-    diffusion.cpp
-    convection_diffusion.cpp
+    #linear_advection.cpp
+    #diffusion.cpp
+    #convection_diffusion.cpp
     physics.cpp
     )
 
diff --git a/src/physics/physics.cpp b/src/physics/physics.cpp
index 21097a6df..b8eb94878 100644
--- a/src/physics/physics.cpp
+++ b/src/physics/physics.cpp
@@ -35,6 +35,20 @@ namespace PHiLiP
     template <int dim, int nstate, typename real>
     Physics<dim,nstate,real>::~Physics() {}
 
+    template <int dim, int nstate, typename real>
+    void Physics<dim,nstate,real>
+    ::dissipative_flux_A_gradu (
+        const real scaling,
+        const std::array<real,nstate> &solution,
+        const std::array<Tensor<1,dim,real>,nstate> &solution_gradient,
+        std::array<Tensor<1,dim,real>,nstate> &dissipative_flux) const
+    {
+        const std::array<Tensor<1,dim,real>,nstate> dissipation = apply_diffusion_matrix(solution, solution_gradient);
+        for (int s=0; s<nstate; s++) {
+            dissipative_flux[s] = -scaling*dissipation[s];
+        }
+    }
+
     // Common manufactured solution for advection, diffusion, convection-diffusion
     template <int dim, int nstate, typename real>
     void Physics<dim,nstate,real>
@@ -51,6 +65,27 @@ namespace PHiLiP
         if (dim==3) uexact[ISTATE] = sin(a*pos[0]+d)*sin(b*pos[1]+e)*sin(c*pos[2]+f);
         solution = uexact;
     }
+    template <int dim, int nstate, typename real>
+    void Physics<dim,nstate,real>
+    ::manufactured_gradient (const Point<dim,double> &pos, std::array<Tensor<1,dim,real>,nstate> &solution_gradient) const
+    {
+        std::array<real,nstate> uexact;
+        
+        using phys = Physics<dim,nstate,real>;
+        const double a = phys::freq_x, b = phys::freq_y, c = phys::freq_z;
+        const double d = phys::offs_x, e = phys::offs_y, f = phys::offs_z;
+        const int ISTATE = 0;
+        if (dim==1) {
+            solution_gradient[ISTATE][0] = a*cos(a*pos[0]+d);
+        } else if (dim==2) {
+            solution_gradient[ISTATE][0] = a*cos(a*pos[0]+d)*sin(b*pos[1]+e);
+            solution_gradient[ISTATE][1] = b*sin(a*pos[0]+d)*cos(b*pos[1]+e);
+        } else if (dim==3) {
+            solution_gradient[ISTATE][0] = a*cos(a*pos[0]+d)*sin(b*pos[1]+e)*sin(c*pos[2]+f);
+            solution_gradient[ISTATE][1] = b*sin(a*pos[0]+d)*cos(b*pos[1]+e)*sin(c*pos[2]+f);
+            solution_gradient[ISTATE][2] = c*sin(a*pos[0]+d)*sin(b*pos[1]+e)*cos(c*pos[2]+f);
+        }
+    }
 
     template <int dim, int nstate, typename real>
     double Physics<dim,nstate,real>
@@ -126,5 +161,282 @@ namespace PHiLiP
     template class PhysicsFactory<PHILIP_DIM, 1, Sacado::Fad::DFad<double> >;
     template class PhysicsFactory<PHILIP_DIM, 1, double>;
 
+    // Linear advection functions
+    template <int dim, int nstate, typename real>
+    Tensor<1,dim,real> LinearAdvection<dim,nstate,real>
+    ::advection_speed () const
+    {
+        Tensor<1,dim,real> advection_speed;
+
+        if(dim >= 1) advection_speed[0] = this->velo_x;
+        if(dim >= 2) advection_speed[1] = this->velo_y;
+        if(dim >= 3) advection_speed[2] = this->velo_z;
+
+        return advection_speed;
+    }
+
+    template <int dim, int nstate, typename real>
+    std::array<real,nstate> LinearAdvection<dim,nstate,real>
+    ::convective_eigenvalues (
+        const std::array<real,nstate> &/*solution*/,
+        const Tensor<1,dim,real> &normal) const
+    {
+        std::array<real,nstate> eig;
+        const Tensor<1,dim,real> advection_speed = this->advection_speed();
+        for (int i=0; i<nstate; i++) {
+            eig[i] = advection_speed*normal;
+        }
+        return eig;
+    }
+
+    template <int dim, int nstate, typename real>
+    void LinearAdvection<dim,nstate,real>
+    ::convective_flux (
+        const std::array<real,nstate> &solution,
+        std::array<Tensor<1,dim,real>,nstate> &conv_flux) const
+    {
+        // Assert conv_flux dimensions
+        const Tensor<1,dim,real> velocity_field = this->advection_speed();
+        for (int i=0; i<nstate; ++i) {
+            conv_flux[i] = velocity_field * solution[i];
+        }
+    }
+
+    template <int dim, int nstate, typename real>
+    std::array<Tensor<1,dim,real>,nstate> LinearAdvection<dim,nstate,real>
+    ::apply_diffusion_matrix(
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_grad) const
+    {
+		std::array<Tensor<1,dim,real>,nstate> zero; // deal.II tensors are initialized with zeros
+        return zero;
+    }
+
+    template <int dim, int nstate, typename real>
+    void LinearAdvection<dim,nstate,real>
+    ::dissipative_flux (
+        const std::array<real,nstate> &solution,
+        const std::array<Tensor<1,dim,real>,nstate> &solution_gradient,
+        std::array<Tensor<1,dim,real>,nstate> &diss_flux) const
+    {
+        // No dissipation
+        const double diff_coeff = this->diff_coeff;
+        for (int i=0; i<nstate; i++) {
+            diss_flux[i] = 0;
+        }
+    }
+
+    template <int dim, int nstate, typename real>
+    void LinearAdvection<dim,nstate,real>
+    ::source_term (
+        const Point<dim,double> &pos,
+        const std::array<real,nstate> &/*solution*/,
+        std::array<real,nstate> &source) const
+    {
+        const Tensor<1,dim,real> vel = this->advection_speed();
+        using phys = Physics<dim,nstate,real>;
+        const double a = phys::freq_x, b = phys::freq_y, c = phys::freq_z;
+        const double d = phys::offs_x, e = phys::offs_y, f = phys::offs_z;
+        const int ISTATE = 0;
+        if (dim==1) {
+            const real x = pos[0];
+            source[ISTATE] = vel[0]*a*cos(a*x+d);
+        } else if (dim==2) {
+            const real x = pos[0], y = pos[1];
+            source[ISTATE] = vel[0]*a*cos(a*x+d)*sin(b*y+e) +
+                             vel[1]*b*sin(a*x+d)*cos(b*y+e);
+        } else if (dim==3) {
+            const real x = pos[0], y = pos[1], z = pos[2];
+            source[ISTATE] =  vel[0]*a*cos(a*x+d)*sin(b*y+e)*sin(c*z+f) +
+                              vel[1]*b*sin(a*x+d)*cos(b*y+e)*sin(c*z+f) +
+                              vel[2]*c*sin(a*x+d)*sin(b*y+e)*cos(c*z+f);
+        }
+    }
+
+    template class LinearAdvection < PHILIP_DIM, 1, Sacado::Fad::DFad<double>  >;
+    template class LinearAdvection < PHILIP_DIM, 1, double >;
+
+    template <int dim, int nstate, typename real>
+    void Diffusion<dim, nstate, real>
+    ::convective_flux (
+        const std::array<real,nstate> &/*solution*/,
+        std::array<Tensor<1,dim,real>,nstate> &/*conv_flux*/) const
+    { }
+
+    template <int dim, int nstate, typename real>
+    std::array<real, nstate> Diffusion<dim, nstate, real>
+    ::convective_eigenvalues(
+        const std::array<real,nstate> &/*solution*/,
+        const Tensor<1,dim,real> &/*normal*/) const
+    {
+        std::array<real,nstate> eig;
+        for (int i=0; i<nstate; i++) {
+            eig[i] = 0;
+        }
+        return eig;
+    }
+
+    template <int dim, int nstate, typename real>
+    std::array<Tensor<1,dim,real>,nstate> Diffusion<dim,nstate,real>
+    ::apply_diffusion_matrix(
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_grad) const
+    {
+        // deal.II tensors are initialized with zeros
+        std::array<Tensor<1,dim,real>,nstate> diffusion;
+        for (int d=0; d<dim; d++) {
+            diffusion[0][d] = 1.0*solution_grad[0][d];
+        }
+        return diffusion;
+    }
+
+    template <int dim, int nstate, typename real>
+    void Diffusion<dim,nstate,real>
+    ::dissipative_flux (
+        const std::array<real,nstate> &solution,
+        const std::array<Tensor<1,dim,real>,nstate> &solution_gradient,
+        std::array<Tensor<1,dim,real>,nstate> &diss_flux) const
+    {
+        const double diff_coeff = this->diff_coeff;
+        for (int i=0; i<nstate; i++) {
+            diss_flux[i] = -diff_coeff*1.0*solution_gradient[i];
+        }
+    }
+
+    template <int dim, int nstate, typename real>
+    void Diffusion<dim,nstate,real>
+    ::source_term (
+        const Point<dim,double> &pos,
+        const std::array<real,nstate> &/*solution*/,
+        std::array<real,nstate> &source) const
+    {
+        using phys = Physics<dim,nstate,real>;
+        const double a = phys::freq_x, b = phys::freq_y, c = phys::freq_z;
+        const double d = phys::offs_x, e = phys::offs_y, f = phys::offs_z;
+        const double diff_coeff = this->diff_coeff;
+        const int ISTATE = 0;
+        if (dim==1) {
+            const real x = pos[0];
+            source[ISTATE] = diff_coeff*a*a*sin(a*x+d);
+        } else if (dim==2) {
+            const real x = pos[0], y = pos[1];
+            source[ISTATE] = diff_coeff*a*a*sin(a*x+d)*sin(b*y+e) +
+                     diff_coeff*b*b*sin(a*x+d)*sin(b*y+e);
+        } else if (dim==3) {
+            const real x = pos[0], y = pos[1], z = pos[2];
+
+            source[ISTATE] =  diff_coeff*a*a*sin(a*x+d)*sin(b*y+e)*sin(c*z+f) +
+                      diff_coeff*b*b*sin(a*x+d)*sin(b*y+e)*sin(c*z+f) +
+                      diff_coeff*c*c*sin(a*x+d)*sin(b*y+e)*sin(c*z+f);
+        }
+    }
+
+    template class Diffusion < PHILIP_DIM, 1, Sacado::Fad::DFad<double>  >;
+    template class Diffusion < PHILIP_DIM, 1, double >;
+
+    template <int dim, int nstate, typename real>
+    void ConvectionDiffusion<dim,nstate,real>
+    ::convective_flux (
+        const std::array<real,nstate> &solution,
+        std::array<Tensor<1,dim,real>,nstate> &conv_flux) const
+    {
+        const Tensor<1,dim,real> velocity_field = this->advection_speed();
+        for (int i=0; i<nstate; ++i) {
+            conv_flux[i] = velocity_field * solution[i];
+        }
+    }
+
+    template <int dim, int nstate, typename real>
+    Tensor<1,dim,real> ConvectionDiffusion<dim,nstate,real>
+    ::advection_speed () const
+    {
+        Tensor<1,dim,real> advection_speed;
+
+        if(dim >= 1) advection_speed[0] = this->velo_x;
+        if(dim >= 2) advection_speed[1] = this->velo_y;
+        if(dim >= 3) advection_speed[2] = this->velo_z;
+        return advection_speed;
+    }
+
+    template <int dim, int nstate, typename real>
+    std::array<real,nstate> ConvectionDiffusion<dim,nstate,real>
+    ::convective_eigenvalues (
+        const std::array<real,nstate> &/*solution*/,
+        const Tensor<1,dim,real> &normal) const
+    {
+        std::array<real,nstate> eig;
+        const Tensor<1,dim,real> advection_speed = this->advection_speed();
+        for (int i=0; i<nstate; i++) {
+            eig[i] = advection_speed*normal;
+        }
+        return eig;
+    }
+
+    template <int dim, int nstate, typename real>
+    std::array<Tensor<1,dim,real>,nstate> ConvectionDiffusion<dim,nstate,real>
+    ::apply_diffusion_matrix(
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_grad) const
+    {
+        // deal.II tensors are initialized with zeros
+        std::array<Tensor<1,dim,real>,nstate> diffusion;
+        for (int d=0; d<dim; d++) {
+            diffusion[0][d] = 1.0*solution_grad[0][d];
+        }
+        return diffusion;
+    }
+
+    template <int dim, int nstate, typename real>
+    void ConvectionDiffusion<dim,nstate,real>
+    ::dissipative_flux (
+        const std::array<real,nstate> &solution,
+        const std::array<Tensor<1,dim,real>,nstate> &solution_gradient,
+        std::array<Tensor<1,dim,real>,nstate> &diss_flux) const
+    {
+        const double diff_coeff = this->diff_coeff;
+        for (int i=0; i<nstate; i++) {
+            diss_flux[i] = -diff_coeff*1.0*solution_gradient[i];
+        }
+    }
+
+    template <int dim, int nstate, typename real>
+    void ConvectionDiffusion<dim,nstate,real>
+    ::source_term (
+        const Point<dim,double> &pos,
+        const std::array<real,nstate> &/*solution*/,
+        std::array<real,nstate> &source) const
+    {
+        const Tensor<1,dim,real> velocity_field = this->advection_speed();
+        using phys = Physics<dim,nstate,real>;
+        const double a = phys::freq_x, b = phys::freq_y, c = phys::freq_z;
+        const double d = phys::offs_x, e = phys::offs_y, f = phys::offs_z;
+
+        const double diff_coeff = this->diff_coeff;
+        const int ISTATE = 0;
+        if (dim==1) {
+            const real x = pos[0];
+            source[ISTATE] = velocity_field[0]*a*cos(a*x+d) +
+                     diff_coeff*a*a*sin(a*x+d);
+        } else if (dim==2) {
+            const real x = pos[0], y = pos[1];
+            source[ISTATE] = velocity_field[0]*a*cos(a*x+d)*sin(b*y+e) +
+                     velocity_field[1]*b*sin(a*x+d)*cos(b*y+e) +
+                     diff_coeff*a*a*sin(a*x+d)*sin(b*y+e) +
+                     diff_coeff*b*b*sin(a*x+d)*sin(b*y+e);
+        } else if (dim==3) {
+            const real x = pos[0], y = pos[1], z = pos[2];
+            source[ISTATE] =   velocity_field[0]*a*cos(a*x+d)*sin(b*y+e)*sin(c*z+f) +
+                       velocity_field[1]*b*sin(a*x+d)*cos(b*y+e)*sin(c*z+f) +
+                       velocity_field[2]*c*sin(a*x+d)*sin(b*y+e)*cos(c*z+f) +
+                       diff_coeff*a*a*sin(a*x+d)*sin(b*y+e)*sin(c*z+f) +
+                       diff_coeff*b*b*sin(a*x+d)*sin(b*y+e)*sin(c*z+f) +
+                       diff_coeff*c*c*sin(a*x+d)*sin(b*y+e)*sin(c*z+f);
+        }
+    }
+    // Instantiate
+    template class ConvectionDiffusion < PHILIP_DIM, 1, double >;
+    template class ConvectionDiffusion < PHILIP_DIM, 1, Sacado::Fad::DFad<double>  >;
+
+
 } // end of PHiLiP namespace
 
diff --git a/src/physics/physics.h b/src/physics/physics.h
index 1d7ff2bd0..ab0755cc4 100644
--- a/src/physics/physics.h
+++ b/src/physics/physics.h
@@ -32,6 +32,9 @@ namespace PHiLiP
         virtual void manufactured_solution (
             const Point<dim,double> &pos,
             std::array<real,nstate> &solution) const;
+        virtual void manufactured_gradient (
+            const Point<dim,double> &pos,
+            std::array<Tensor<1,dim,real>,nstate> &solution_gradient) const;
 
         /// Returns the integral of the manufactured solution over the hypercube [0,1]
         ///
@@ -51,8 +54,17 @@ namespace PHiLiP
             const Tensor<1,dim,real> &/*normal*/) const = 0;
 
         // Evaluate the diffusion matrix $ A $ such that $F_v = A \nabla u$
-        virtual std::array<Tensor<2,dim,real>,nstate> diffusion_matrix (
-            const std::array<real,nstate> &solution) const = 0;
+        virtual std::array<Tensor<1,dim,real>,nstate> apply_diffusion_matrix (
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_grad) const = 0;
+
+        // Dissipative fluxes that will be differentiated once in space
+        // Evaluates the dissipative flux through the linearization F = A(u)*grad(u)
+        void dissipative_flux_A_gradu (
+            const real scaling,
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_gradient,
+            std::array<Tensor<1,dim,real>,nstate> &diss_flux) const;
 
         // Dissipative fluxes that will be differentiated once in space
         virtual void dissipative_flux (
@@ -92,7 +104,7 @@ namespace PHiLiP
         //const double offs_x = 1, offs_y = 1.2, offs_z = 1.5;
 
         const double pi = atan(1)*4.0;
-        const double freq_x = 1.59/dim, freq_y = 2*1.81/dim,    freq_z = 3*1.76/dim;
+        const double freq_x = 0.59/dim, freq_y = 2*0.81/dim,    freq_z = 3*0.76/dim;
         const double offs_x = 1,        offs_y = 1.2,           offs_z = 1.5;
         const double velo_x = exp(1)/2, velo_y =-pi/4.0,        velo_z = sqrt(2);
         //const double velo_x = 1.0, velo_y =-pi/4.0,        velo_z = sqrt(2);
@@ -129,8 +141,9 @@ namespace PHiLiP
             const Tensor<1,dim,real> &normal) const;
 
         // Diffusion matrix: 0
-        std::array<Tensor<2,dim,real>,nstate> diffusion_matrix (
-            const std::array<real,nstate> &solution) const;
+        std::array<Tensor<1,dim,real>,nstate> apply_diffusion_matrix (
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_grad) const;
 
         // Dissipative flux: 0
         void dissipative_flux (
@@ -171,8 +184,9 @@ namespace PHiLiP
             const Tensor<1,dim,real> &/*normal*/) const;
 
         // Diffusion matrix is identity
-        std::array<Tensor<2,dim,real>,nstate> diffusion_matrix (
-            const std::array<real,nstate> &solution) const;
+        std::array<Tensor<1,dim,real>,nstate> apply_diffusion_matrix (
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_grad) const;
 
         // Dissipative flux: u
         void dissipative_flux (
@@ -208,8 +222,9 @@ namespace PHiLiP
             const Tensor<1,dim,real> &/*normal*/) const;
 
         // Diffusion matrix is identity
-        std::array<Tensor<2,dim,real>,nstate> diffusion_matrix (
-            const std::array<real,nstate> &solution) const;
+        std::array<Tensor<1,dim,real>,nstate> apply_diffusion_matrix (
+            const std::array<real,nstate> &solution,
+            const std::array<Tensor<1,dim,real>,nstate> &solution_grad) const;
 
         // Dissipative flux: u
         void dissipative_flux (
diff --git a/tests/convection_diffusion_implicit/1d_convection_diffusion_implicit.prm b/tests/convection_diffusion_implicit/1d_convection_diffusion_implicit.prm
index 4a761af40..45df15016 100644
--- a/tests/convection_diffusion_implicit/1d_convection_diffusion_implicit.prm
+++ b/tests/convection_diffusion_implicit/1d_convection_diffusion_implicit.prm
@@ -4,15 +4,18 @@
 set dimension = 1
 
 # The PDE we want to solve. Choices are
-# <convection_diffusion|diffusion|convection_diffusion>.
+# <diffusion|diffusion|convection_diffusion>.
 set pde_type  = convection_diffusion
 
 subsection ODE solver
+
+  set ode_output                          = verbose
+
   # Maximum nonlinear solver iterations
   set nonlinear_max_iterations            = 500000
 
   # Nonlinear solver residual tolerance
-  set nonlinear_steady_residual_tolerance = 1e-13
+  set nonlinear_steady_residual_tolerance = 1e-12
 
   # Print every print_iteration_modulo iterations of the nonlinear solver
   set print_iteration_modulo              = 1
@@ -26,15 +29,15 @@ subsection manufactured solution convergence study
   set degree_end        = 3
 
   # Starting degree for convergence study
-  set degree_start      = 0
+  set degree_start      = 1
 
   # Multiplier on grid size. nth-grid will be of size
   # (initial_grid^grid_progression)^dim
   set grid_progression  = 1.5
 
   # Initial grid of size (initial_grid_size)^dim
-  set initial_grid_size = 2
+  set initial_grid_size = 3
 
   # Number of grids in grid study
-  set number_of_grids   = 4
+  set number_of_grids   = 6
 end
diff --git a/tests/convection_diffusion_implicit/2d_convection_diffusion_implicit.prm b/tests/convection_diffusion_implicit/2d_convection_diffusion_implicit.prm
index 4f6370bec..7a926b2ac 100644
--- a/tests/convection_diffusion_implicit/2d_convection_diffusion_implicit.prm
+++ b/tests/convection_diffusion_implicit/2d_convection_diffusion_implicit.prm
@@ -5,7 +5,7 @@ set dimension = 2
 
 # The PDE we want to solve. Choices are
 # <advection|diffusion|convection_diffusion>.
-set pde_type  = advection
+set pde_type  = convection_diffusion
 
 subsection ODE solver
   # Maximum nonlinear solver iterations
@@ -26,7 +26,7 @@ subsection manufactured solution convergence study
   set degree_end        = 3
 
   # Starting degree for convergence study
-  set degree_start      = 0
+  set degree_start      = 1
 
   # Multiplier on grid size. nth-grid will be of size
   # (initial_grid^grid_progression)^dim
@@ -36,6 +36,6 @@ subsection manufactured solution convergence study
   set initial_grid_size = 2
 
   # Number of grids in grid study
-  set number_of_grids   = 4
+  set number_of_grids   = 5
 end
 
diff --git a/tests/convection_diffusion_implicit/3d_convection_diffusion_implicit.prm b/tests/convection_diffusion_implicit/3d_convection_diffusion_implicit.prm
index 95026a42f..46a117c34 100644
--- a/tests/convection_diffusion_implicit/3d_convection_diffusion_implicit.prm
+++ b/tests/convection_diffusion_implicit/3d_convection_diffusion_implicit.prm
@@ -5,7 +5,7 @@ set dimension = 3
 
 # The PDE we want to solve. Choices are
 # <advection|diffusion|convection_diffusion>.
-set pde_type  = advection
+set pde_type  = convection_diffusion
 
 subsection ODE solver
   # Maximum nonlinear solver iterations
@@ -26,7 +26,7 @@ subsection manufactured solution convergence study
   set degree_end        = 3
 
   # Starting degree for convergence study
-  set degree_start      = 0
+  set degree_start      = 1
 
   # Multiplier on grid size. nth-grid will be of size
   # (initial_grid^grid_progression)^dim
diff --git a/tests/diffusion_implicit/2d_diffusion_implicit.prm b/tests/diffusion_implicit/2d_diffusion_implicit.prm
index 68c515662..b11f2ca19 100644
--- a/tests/diffusion_implicit/2d_diffusion_implicit.prm
+++ b/tests/diffusion_implicit/2d_diffusion_implicit.prm
@@ -26,7 +26,7 @@ subsection manufactured solution convergence study
   set degree_end        = 3
 
   # Starting degree for convergence study
-  set degree_start      = 2
+  set degree_start      = 1
 
   # Multiplier on grid size. nth-grid will be of size
   # (initial_grid^grid_progression)^dim
@@ -36,6 +36,6 @@ subsection manufactured solution convergence study
   set initial_grid_size = 2
 
   # Number of grids in grid study
-  set number_of_grids   = 4
+  set number_of_grids   = 5
 end
 
diff --git a/tests/numerical_flux/numerical_flux_conservation.cpp b/tests/numerical_flux/numerical_flux_conservation.cpp
index b54aa94d6..5de493ac0 100644
--- a/tests/numerical_flux/numerical_flux_conservation.cpp
+++ b/tests/numerical_flux/numerical_flux_conservation.cpp
@@ -1,3 +1,6 @@
+#include <deal.II/base/tensor.h>
+
+#include "parameters.h"
 #include "parameters.h"
 #include "physics/physics.h"
 #include "numerical_flux/numerical_flux.h"
@@ -6,65 +9,148 @@ using PDEType  = Parameters::AllParameters::PartialDifferentialEquation;
 using ConvType = Parameters::AllParameters::ConvectiveNumericalFlux;
 using DissType = Parameters::AllParameters::DissipativeNumericalFlux;
 
+const double TOLERANCE = 1E-12;
+
 template<int dim, int nstate>
-int test_numerical_flux_conservation (const PDEType pde_type, const ConvType conv_type)
+void set_solution (
+    std::array<double, nstate>               &soln_int, 
+    std::array<double, nstate>               &soln_ext,
+    std::array<dealii::Tensor<1,dim,double>, nstate> &soln_grad_int,
+    std::array<dealii::Tensor<1,dim,double>, nstate> &soln_grad_ext,
+    dealii::Tensor<1,dim,double>                     &normal_int)
 {
-    const double TOLERANCE = 1E-12;
-    std::cout << std::setprecision(std::numeric_limits<double>::digits10 + 1);
-    std::cout << std::scientific;
-
-
-    using namespace PHiLiP;
-    Physics<dim, nstate, double> *pde_physics = PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
-
-    NumericalFluxConvective<dim, nstate, double> *conv_num_flux = 
-        NumericalFluxFactory<dim, nstate, double>
-        ::create_convective_numerical_flux (conv_type, pde_physics);
-
-    Tensor<1,dim,double> normal_int;
     for (int d=0; d<dim; d++) {
         normal_int[d] = 1;
     }
 
-    std::array<double, nstate> soln_int, soln_ext;
     for (int s=0; s<nstate; s++) {
         soln_int[s] = exp(s+exp(1));
         soln_ext[s] = sin(s-atan(1));
-    }
 
-    std::array<double, nstate> conv_num_flux_dot_n_1 = conv_num_flux->evaluate_flux(soln_int, soln_ext, normal_int);
-    std::array<double, nstate> conv_num_flux_dot_n_2 = conv_num_flux->evaluate_flux(soln_ext, soln_int, -normal_int);
+        for (int d=0; d<dim; d++) {
+            soln_int[s] = d*exp(s+d*exp(1));
+            soln_ext[s] = sin(d*s-atan(1))/(d+0.1);
+        }
+    }
+}
 
-    std::cout << "PDE type " << pde_type << std::endl;
+template<int dim, int nstate>
+void compare_array ( const std::array<double, nstate> &array1, const std::array<double, nstate> &array2, double scale2)
+{
     for (int s=0; s<nstate; s++) {
-        // Flux should be equal and opposite, therefore f1 + f2 = 0
-        const double diff = std::abs(conv_num_flux_dot_n_1[s] + conv_num_flux_dot_n_2[s]);
+        const double diff = std::abs(array1[s] - scale2*array2[s]);
         std::cout
             << "State " << s << " out of " << nstate
             << std::endl
-            << "Num flux 1 = " << conv_num_flux_dot_n_1[s]
+            << "Array 1 = " << array1[s]
             << std::endl
-            << "Num flux 2 = " << conv_num_flux_dot_n_2[s]
+            << "Array 2 = " << array2[s]
             << std::endl
             << "Difference = " << diff
             << std::endl;
         assert(diff < TOLERANCE);
     }
+    std::cout << std::endl
+              << std::endl
+              << std::endl;
+}
 
-    delete conv_num_flux;
+template<int dim, int nstate>
+int test_dissipative_numerical_flux_conservation (const PDEType pde_type, const DissType diss_type)
+{
+    std::cout << std::setprecision(std::numeric_limits<double>::digits10 + 1);
+    std::cout << std::scientific;
+
+
+    using namespace PHiLiP;
+    Physics<dim, nstate, double> *pde_physics = PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
+
+    Tensor<1,dim,double> normal_int;
+    std::array<double, nstate> soln_int, soln_ext;
+    std::array<Tensor<1,dim,double>, nstate> soln_grad_int, soln_grad_ext;
+    set_solution<dim, nstate> (soln_int, soln_ext, soln_grad_int, soln_grad_ext, normal_int);
+
+    NumericalFluxDissipative<dim, nstate, double> *diss_num_flux = 
+        NumericalFluxFactory<dim, nstate, double>
+        ::create_dissipative_numerical_flux (diss_type, pde_physics);
+    std::array<double, nstate> diss_num_flux_dot_n_1 = diss_num_flux->evaluate_solution_flux(soln_int, soln_ext, normal_int);
+    std::array<double, nstate> diss_num_flux_dot_n_2 = diss_num_flux->evaluate_solution_flux(soln_ext, soln_int, -normal_int);
+
+    double penalty = 100;
+    std::array<double, nstate> diss_auxi_num_flux_dot_n_1 = diss_num_flux->evaluate_auxiliary_flux(
+                 soln_int, soln_ext,
+                 soln_grad_int, soln_grad_ext,
+                 normal_int, penalty);
+
+    std::array<double, nstate> diss_auxi_num_flux_dot_n_2 = diss_num_flux->evaluate_auxiliary_flux(
+                 soln_ext, soln_int,
+                 soln_grad_ext, soln_grad_int,
+                 -normal_int, penalty);
+ 
+
+    std::cout << "Dissipative solution flux conservation (should be equal since not dotted with normal)" << std::endl;
+    compare_array<dim,nstate> (diss_num_flux_dot_n_1, diss_num_flux_dot_n_2, 1.0);
+
+    std::cout << "Dissipative auxiliary flux conservation (should be equal and opposite since it is dotted with normal)" << std::endl;
+    compare_array<dim,nstate> (diss_auxi_num_flux_dot_n_1, diss_auxi_num_flux_dot_n_2, -1.0);
+
+    delete diss_num_flux;
     delete pde_physics;
 
     return 0;
 }
 
 template<int dim, int nstate>
-int test_numerical_flux_consistency (const PDEType pde_type, const ConvType conv_type)
+int test_dissipative_numerical_flux_consistency (const PDEType pde_type, const DissType diss_type)
 {
-    const double TOLERANCE = 1E-12;
-    std::cout << std::setprecision(std::numeric_limits<double>::digits10 + 1);
-    std::cout << std::scientific;
+    using namespace PHiLiP;
+    Physics<dim, nstate, double> *pde_physics = PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
+
+    NumericalFluxDissipative<dim, nstate, double> *diss_num_flux = 
+        NumericalFluxFactory<dim, nstate, double>
+        ::create_dissipative_numerical_flux (diss_type, pde_physics);
+
+    Tensor<1,dim,double> normal_int;
+    std::array<double, nstate> soln_int, soln_ext;
+    std::array<Tensor<1,dim,double>, nstate> soln_grad_int, soln_grad_ext;
+    set_solution<dim,nstate> (soln_int, soln_ext, soln_grad_int, soln_grad_ext, normal_int);
+
+    // Copy state to ext cell
+    soln_int = soln_ext;
+    soln_grad_int = soln_grad_ext;
 
+    // Evaluate numerical fluxes
+    std::array<double, nstate> diss_soln_num_flux_dot_n = diss_num_flux->evaluate_flux(soln_int, soln_ext, normal_int);
+    double penalty = 100;
+    std::array<double, nstate> diss_auxi_num_flux_dot_n = diss_num_flux->evaluate_auxiliary_flux(
+                 soln_int, soln_ext,
+                 soln_grad_int, soln_grad_ext,
+                 normal_int, penalty);
 
+
+    std::array<Tensor<1,dim,double>, nstate> diss_phys_flux_int;
+    pde_physics->dissipative_flux (soln_int, diss_phys_flux_int);
+
+    std::array<double, nstate> diss_phys_flux_dot_n;
+    for (int s=0; s<nstate; s++) {
+        diss_phys_flux_dot_n[s] = diss_phys_flux_int[s]*normal_int;
+    }
+
+    std::cout << "Dissipative solution flux should be consistent (equal) with solution, when same values on each side of boundary" << std::endl;
+    compare_array<dim,nstate> (diss_soln_num_flux_dot_n, soln_int, 1.0);
+
+    std::cout << "Dissipative auxiliary flux should be consistent (equal) with gradient, when same values on each side of boundary" << std::endl;
+    compare_array<dim,nstate> (diss_auxi_num_flux_dot_n, diss_phys_flux_int, 1.0);
+
+    delete diss_num_flux;
+    delete pde_physics;
+
+    return 0;
+}
+
+template<int dim, int nstate>
+int test_convective_numerical_flux_conservation (const PDEType pde_type, const ConvType conv_type)
+{
     using namespace PHiLiP;
     Physics<dim, nstate, double> *pde_physics = PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
 
@@ -73,16 +159,40 @@ int test_numerical_flux_consistency (const PDEType pde_type, const ConvType conv
         ::create_convective_numerical_flux (conv_type, pde_physics);
 
     Tensor<1,dim,double> normal_int;
-    for (int d=0; d<dim; d++) {
-        normal_int[d] = 1;
-    }
+    std::array<double, nstate> soln_int, soln_ext;
+    std::array<Tensor<1,dim,double>, nstate> soln_grad_int, soln_grad_ext;
+    set_solution<dim,nstate> (soln_int, soln_ext, soln_grad_int, soln_grad_ext, normal_int);
+
+    std::array<double, nstate> conv_num_flux_dot_n_1 = conv_num_flux->evaluate_flux(soln_int, soln_ext, normal_int);
+    std::array<double, nstate> conv_num_flux_dot_n_2 = conv_num_flux->evaluate_flux(soln_ext, soln_int, -normal_int);
+
+    std::cout << "Convective numerical flux conservation (should be equal and opposite)" << std::endl;
+    compare_array<dim,nstate> (conv_num_flux_dot_n_1, conv_num_flux_dot_n_2, -1.0);
+
+    delete conv_num_flux;
+    delete pde_physics;
+
+    return 0;
+}
 
+template<int dim, int nstate>
+int test_convective_numerical_flux_consistency (const PDEType pde_type, const ConvType conv_type)
+{
+    using namespace PHiLiP;
+    Physics<dim, nstate, double> *pde_physics = PhysicsFactory<dim, nstate, double>::create_Physics(pde_type);
+
+    NumericalFluxConvective<dim, nstate, double> *conv_num_flux = 
+        NumericalFluxFactory<dim, nstate, double>
+        ::create_convective_numerical_flux (conv_type, pde_physics);
+
+    Tensor<1,dim,double> normal_int;
     std::array<double, nstate> soln_int, soln_ext;
-    for (int s=0; s<nstate; s++) {
-        soln_int[s] = sin(s+exp(1));
-        soln_ext[s] = soln_int[s];
-    }
+    std::array<Tensor<1,dim,double>, nstate> soln_grad_int, soln_grad_ext;
+    set_solution<dim,nstate> (soln_int, soln_ext, soln_grad_int, soln_grad_ext, normal_int);
 
+    // Consistent numerical flux should be equal to physical flux when both states are equal
+    // Therefore, f1 - f2 = 0
+    soln_int = soln_ext;
     std::array<double, nstate> conv_num_flux_dot_n = conv_num_flux->evaluate_flux(soln_int, soln_ext, normal_int);
 
     std::array<Tensor<1,dim,double>, nstate> conv_phys_flux_int;
@@ -93,26 +203,8 @@ int test_numerical_flux_consistency (const PDEType pde_type, const ConvType conv
         conv_phys_flux_dot_n[s] = conv_phys_flux_int[s]*normal_int;
     }
 
-
-    std::cout << "PDE type " << pde_type << std::endl;
-    for (int s=0; s<nstate; s++) {
-        // Consistent numerical flux should be equal to physical flux when both states are equal
-        // Therefore, f1 - f2 = 0
-        const double diff = std::abs(conv_num_flux_dot_n[s] - conv_phys_flux_dot_n[s]);
-        std::cout
-            << "State " << s << " out of " << nstate
-            << std::endl
-            << "Num flux 1 = " << conv_num_flux_dot_n[s]
-            << std::endl
-            << "Num flux 2 = " << conv_phys_flux_dot_n[s]
-            << std::endl
-            << "Difference = " << diff
-            << std::endl;
-        assert(diff < TOLERANCE);
-    }
-    std::cout << std::endl
-              << std::endl
-              << std::endl;
+    std::cout << "Convective numerical flux should be consistent with physical flux (equal), when same values on each side of boundary" << std::endl;
+    compare_array<dim,nstate> (conv_num_flux_dot_n, conv_phys_flux_dot_n, 1.0);
 
     delete conv_num_flux;
     delete pde_physics;
@@ -129,6 +221,7 @@ int main (int argc, char *argv[])
         PDEType::diffusion,
         PDEType::convection_diffusion
     };
+
     std::vector<ConvType> conv_type {
         ConvType::lax_friedrichs
     };
@@ -139,8 +232,10 @@ int main (int argc, char *argv[])
     int success = 0;
     for (auto pde = pde_type.begin(); pde != pde_type.end() && success == 0; pde++) {
         for (auto conv = conv_type.begin(); conv != conv_type.end() && success == 0; conv++) {
-            success = test_numerical_flux_conservation<PHILIP_DIM,1> (*pde, *conv);
-            success = test_numerical_flux_consistency<PHILIP_DIM,1> (*pde, *conv);
+            success = test_convective_numerical_flux_conservation<PHILIP_DIM,1> (*pde, *conv);
+        }
+        for (auto diss = diss_type.begin(); diss != diss_type.end() && success == 0; diss++) {
+            success = test_dissipative_numerical_flux_conservation<PHILIP_DIM,1> (*pde, *diss);
         }
     }
     return success;