/* Distributed under the Apache License, Version 2.0.
   See accompanying NOTICE file for details.*/

#define _USE_MATH_DEFINES
#include "LBMVerification.h"
#include "FileUtils.h"
#include <iostream>
#include <cmath>
#include <limits>
#include <fstream>
#ifdef VTK_EXTENSION
#include "VTKUtils.h"
#endif
#include "RunUtils.h"

void LBM::Verification::RunComparison(std::string const& baseline_root)
{
  auto& status = m_verification_report.m_run_status[m_run_name];

  ///////////////////////////////
  // Run The Analytic Solution //
  ///////////////////////////////

  if (!AnalyticSolution(m_analytic_run))
  {
    status.errors.push_back("Unable to run Analytic Solution.");
    return;
  }
#ifdef VTK_EXTENSION
  Info("Writing Analytic VTR file...");
  LBM::VTK::WriteVTR(m_out_dir + m_run_name + "_analytic.vtr", m_analytic_run);
#endif

  ////////////////////////////////////////
  // Compute Analytic Velocity Profiles //
  ////////////////////////////////////////

  // Load baselines
  bool have_analytic_velocity_baseline = true;
  LBM::VelocityProfiles computed_analytic_velocity_profiles;
  LBM::VelocityProfiles baseline_analytic_velocity_profiles;
  have_analytic_velocity_baseline &= LBM::ReadCSV(baseline_root + m_run_name + "_analytic_y_axis_velocity_profiles_baseline.csv", baseline_analytic_velocity_profiles.y);
  have_analytic_velocity_baseline &= LBM::ReadCSV(baseline_root + m_run_name + "_analytic_z_axis_velocity_profiles_baseline.csv", baseline_analytic_velocity_profiles.z);
  GetVelocityProfiles(m_analytic_run, computed_analytic_velocity_profiles, 40);
  // Compare and write out results
  LBM::CSV analytic_y_velocity_error;
  LBM::CSV analytic_z_velocity_error;
  if (have_analytic_velocity_baseline)
  {
    m_verification_report.CompareProfiles(VerificationReport::ComparisonTypes::AnalyticVelocity,
      m_run_name, "Computed Analytic vs. Baseline Analytic Y Velocity Profiles", baseline_analytic_velocity_profiles.y, computed_analytic_velocity_profiles.y, analytic_y_velocity_error);
    m_verification_report.CompareProfiles(VerificationReport::ComparisonTypes::AnalyticVelocity,
      m_run_name, "Computed Analytic vs. Baseline Analytic Z Velocity Profiles", baseline_analytic_velocity_profiles.z, computed_analytic_velocity_profiles.z, analytic_z_velocity_error);
  }
  else
  {
    m_verification_report.m_run_profile_comparisons[m_run_name]["MISSING Analytic LBM vs. Baseline Analytic Y Velocity Profiles"];
    m_verification_report.m_run_profile_comparisons[m_run_name]["MISSING Analytic LBM vs. Baseline Analytic Z Velocity Profiles"];
    status.errors.push_back("Analytic velocity baseline not found");
    status.analytic_velocity = false;
  }
  // Write out csv as the new baseline, if needed
  Info("Writing Analytic Baseline CSV files...");
  LBM::WriteCSV(m_out_dir + m_run_name + "_analytic_y_axis_velocity_profiles_baseline.csv", computed_analytic_velocity_profiles.y);
  LBM::WriteCSV(m_out_dir + m_run_name + "_analytic_z_axis_velocity_profiles_baseline.csv", computed_analytic_velocity_profiles.z);

  ///////////////////////////////////////////////////////
  // Compare Pressure, Temperature, Wall Stress values //
  ///////////////////////////////////////////////////////

  LBM::CSV analytic_pressure_temp_wall_stress_baseline;
  LBM::CSV computed_analytic_pressure_temperature_wall_stress_values;
  LBM::CSV computed_analytic_pressure_temperature_wall_stress_error;
  // Load baselines
  bool have_analytic_pressure_temp_wall_stress_baseline = LBM::ReadCSV(baseline_root + m_run_name + "_analytic_pressure_temperature_wall_stress_baseline.csv", analytic_pressure_temp_wall_stress_baseline);
  // Compare and write out results
  GetPressureTemperatureWallStressValues(m_analytic_run, computed_analytic_pressure_temperature_wall_stress_values, true);
  if (have_analytic_pressure_temp_wall_stress_baseline)
  {
    m_verification_report.CompareValues(m_run_name, 
                                       "Computed Analytic vs. Baseline Analytic Pressure, Temperature, Wall Stress Values",
                                        analytic_pressure_temp_wall_stress_baseline,
                                        computed_analytic_pressure_temperature_wall_stress_values,
                                        computed_analytic_pressure_temperature_wall_stress_error);
  }
  else
  {
    m_verification_report.m_run_profile_comparisons[m_run_name]["MISSING Analytic LBM vs. Baseline Analytic Pressure, Temperature, Wall Stress"];
    status.errors.push_back("Analytic pressure, temperature, wall stress baseline not found");
    status.analytic_temperature = false;
  }
  // Write out csv as the new baseline, if needed
  Info("Writing Analytic Pressure, Temperature, Wall Stress Baseline CSV files...");
  LBM::WriteCSV(m_out_dir + m_run_name + "_analytic_pressure_temperature_wall_stress_baseline.csv", computed_analytic_pressure_temperature_wall_stress_values);

  //////////////////
  // Run the LBM ///
  //////////////////

  Info("Running LBM Solution");
  if (!m_lbm_run.Compute())
  {
    status.errors.push_back("Error running LBM.");
    return;
  }
  status.num_nodes = m_lbm_run.out.num_cells;
  status.lbm_compute_time_s = m_lbm_run.out.compute_time_s;
  status.num_preconditioner_iterations = m_lbm_run.out.num_preconditioner_iterations;
  status.num_fluid_iterations = m_lbm_run.out.num_fluid_iterations;
  status.num_thermal_iterations = m_lbm_run.out.num_thermal_iterations;
#ifdef VTK_EXTENSION
  Info("Writing LBM VTR file...");
  LBM::VTK::WriteVTR(m_out_dir + m_run_name + "_lbm.vtr", m_lbm_run);
#endif

  ///////////////////////////////////
  // Compute LBM Velocity Profiles //
  ///////////////////////////////////

  Info("Writing LBM CSV files...");

  ////////////////////
  // Velocity Plots //
  ////////////////////

  ///////////////////////////////
  // Velocity Error and Timing //
  ///////////////////////////////

   // Preconditions error/timing
  LBM::CSV precondition_error;
  LBM::CSV precondition_timing;
  std::vector<float>  iterations;
  if (!m_lbm_run.out.pre_conditioning_iteration_error.empty())
  {
    for (size_t i = 0; i < m_lbm_run.out.pre_conditioning_iteration_error.size(); i++)
      iterations.push_back(i);
    precondition_error.push_back(std::pair<std::string, std::vector<float>>("iteration", iterations));
    precondition_error.push_back(std::pair<std::string, std::vector<float>>("error", m_lbm_run.out.pre_conditioning_iteration_error));
    LBM::WriteCSV(m_out_dir + m_run_name + "_precondition_iteration_error.csv", precondition_error);
    precondition_timing.push_back(std::pair<std::string, std::vector<float>>("iteration", iterations));
    precondition_timing.push_back(std::pair<std::string, std::vector<float>>("time_ms", m_lbm_run.out.pre_conditioning_iteration_timing_ms));
    LBM::WriteCSV(m_out_dir + m_run_name + "_precondition_iteration_timing.csv", precondition_timing);
  }

  LBM::CSV velocity_error;
  LBM::CSV velocity_timing;
  iterations.clear();
  for (size_t i = 0; i < m_lbm_run.out.fluid_iteration_error.size(); i++)
    iterations.push_back(i);
  velocity_error.push_back(std::pair<std::string, std::vector<float>>("iteration", iterations));
  velocity_error.push_back(std::pair<std::string, std::vector<float>>("error", m_lbm_run.out.fluid_iteration_error));
  LBM::WriteCSV(m_out_dir + m_run_name + "_pressure_flow_stress_iteration_error.csv", velocity_error);
  velocity_timing.push_back(std::pair<std::string, std::vector<float>>("iteration", iterations));
  velocity_timing.push_back(std::pair<std::string, std::vector<float>>("time_ms", m_lbm_run.out.fluid_iteration_timing_ms));
  LBM::WriteCSV(m_out_dir + m_run_name + "_pressure_flow_stress_iteration_timing.csv", velocity_timing);

#ifdef VTK_EXTENSION
  // Plot the iteration errors
  if (!m_lbm_run.out.pre_conditioning_iteration_error.empty())
  {
    std::vector<VTK::CSVPlotSource> precondition_error_source;
    precondition_error_source.push_back(VTK::CSVPlotSource("Precondition Iteration Error", precondition_error, 1.0, 0.0, 0.0, 2.f, 1/*SOLID_LINE*/));
    LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_precondition_", precondition_error_source, true, *GetLogger());
    std::vector<VTK::CSVPlotSource> precondition_timing_source;
    precondition_timing_source.push_back(VTK::CSVPlotSource("Precondition Iteration Time(ms)", precondition_timing, 1.0, 0.0, 0.0, 2.f, 1/*SOLID_LINE*/));
    LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_precondition_", precondition_timing_source, true, *GetLogger());
  }
  // Pressure, Flow, and Stress
  std::vector<VTK::CSVPlotSource> fluid_error_source;
  fluid_error_source.push_back(VTK::CSVPlotSource("Velocity Iteration Error", velocity_error, 1.0, 0.0, 0.0, 2.f, 1/*SOLID_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_fluid_", fluid_error_source, true, *GetLogger());
  std::vector<VTK::CSVPlotSource> fluid_timing_sources;
  fluid_timing_sources.push_back(VTK::CSVPlotSource("Velocity Iteration Time(ms)", velocity_timing, 1.0, 0.0, 0.0, 2.f, 1/*SOLID_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_fluid_", fluid_timing_sources, true, *GetLogger());
#endif

  ////////////////////////////////////////////////
  // Compare Analytic and LBM Velocity Profiles //
  ////////////////////////////////////////////////

  LBM::CSV lbm_y_velocity_error;
  LBM::CSV lbm_z_velocity_error;
  bool have_lbm_velocity_baseline = true;
  LBM::VelocityProfiles baseline_lbm_velocity_profiles;
  LBM::VelocityProfiles computed_lbm_velocity_profiles;
  have_lbm_velocity_baseline &= LBM::ReadCSV(baseline_root + m_run_name + "_lbm_y_axis_velocity_profiles_baseline.csv", baseline_lbm_velocity_profiles.y);
  have_lbm_velocity_baseline &= LBM::ReadCSV(baseline_root + m_run_name + "_lbm_z_axis_velocity_profiles_baseline.csv", baseline_lbm_velocity_profiles.z);
  GetVelocityProfiles(m_lbm_run, computed_lbm_velocity_profiles, 40);
  if (have_lbm_velocity_baseline)
  {
    m_verification_report.CompareProfiles(VerificationReport::ComparisonTypes::VS_BaselineVelocity,
      m_run_name, "Computed LBM vs. Baseline LBM Y Velocity Profiles", baseline_lbm_velocity_profiles.y, computed_lbm_velocity_profiles.y, lbm_y_velocity_error);
    m_verification_report.CompareProfiles(VerificationReport::ComparisonTypes::VS_BaselineVelocity,
      m_run_name, "Computed LBM vs. Baseline LBM Z Velocity Profiles", baseline_lbm_velocity_profiles.z, computed_lbm_velocity_profiles.z, lbm_z_velocity_error);
  }
  else
  {
    m_verification_report.m_run_profile_comparisons[m_run_name]["MISSING Computed LBM vs. Baseline LBM Y Velocity Profiles"];
    m_verification_report.m_run_profile_comparisons[m_run_name]["MISSING Computed LBM vs. Baseline LBM Z Velocity Profiles"];
    status.errors.push_back("LBM velocity baseline not found");
    status.vs_baseline_velocity = false;
  }
  // Write out csv as the new baseline, if needed
  LBM::WriteCSV(m_out_dir + m_run_name + "_lbm_y_axis_velocity_profiles_baseline.csv", computed_lbm_velocity_profiles.y);
  LBM::WriteCSV(m_out_dir + m_run_name + "_lbm_z_axis_velocity_profiles_baseline.csv", computed_lbm_velocity_profiles.z);

  // Join the analytic and LBM CSVs into one to write
  LBM::CSV y_axis_velocity_profiles;
  for (size_t i = 0; i< computed_analytic_velocity_profiles.y.size(); i++)
  {
    if(i==0)
      y_axis_velocity_profiles.push_back(computed_analytic_velocity_profiles.y[i]);
    else
    {
      y_axis_velocity_profiles.push_back(computed_analytic_velocity_profiles.y[i]);
      y_axis_velocity_profiles.back().first = "Analytic_Velocity_" + computed_analytic_velocity_profiles.y[i].first;
      y_axis_velocity_profiles.push_back(computed_lbm_velocity_profiles.y[i]);
      y_axis_velocity_profiles.back().first = "Computed_LBM_Velocity_" + computed_lbm_velocity_profiles.y[i].first;
      if (have_lbm_velocity_baseline)
      {
        y_axis_velocity_profiles.push_back(baseline_lbm_velocity_profiles.y[i]);
        y_axis_velocity_profiles.back().first = "Baseline_LBM_Velocity_" + baseline_lbm_velocity_profiles.y[i].first;
      }
    }
  }
  LBM::CSV z_axis_velocity_profiles;
  for (size_t i = 0; i< computed_analytic_velocity_profiles.z.size(); i++)
  {
    if (i == 0)
      z_axis_velocity_profiles.push_back(computed_analytic_velocity_profiles.z[i]);
    else
    {
      z_axis_velocity_profiles.push_back(computed_analytic_velocity_profiles.z[i]);
      z_axis_velocity_profiles.back().first = "Analytic_Velocity_" + computed_analytic_velocity_profiles.z[i].first;
      z_axis_velocity_profiles.push_back(computed_lbm_velocity_profiles.z[i]);
      z_axis_velocity_profiles.back().first = "Computed_LBM_Velocity_" + computed_lbm_velocity_profiles.z[i].first;
      if (have_lbm_velocity_baseline)
      {
        z_axis_velocity_profiles.push_back(baseline_lbm_velocity_profiles.z[i]);
        z_axis_velocity_profiles.back().first = "Baseline_LBM_Velocity_" + baseline_lbm_velocity_profiles.z[i].first;
      }
    }
  }
  Info("Writing Full CSV files...");
  LBM::WriteCSV(m_out_dir + m_run_name + "_y_axis_velocity_profiles.csv", y_axis_velocity_profiles);
  LBM::WriteCSV(m_out_dir + m_run_name + "_z_axis_velocity_profiles.csv", z_axis_velocity_profiles);
  // Compare and write out results
  LBM::CSV y_velocity_error;
  LBM::CSV z_velocity_error;
  m_verification_report.CompareProfiles(VerificationReport::ComparisonTypes::VS_AnalyticVelocity,
    m_run_name,"Computed LBM vs. Analytic Y Velocity Profiles", computed_analytic_velocity_profiles.y, computed_lbm_velocity_profiles.y, y_velocity_error);
  m_verification_report.CompareProfiles(VerificationReport::ComparisonTypes::VS_AnalyticVelocity,
    m_run_name,"Computed LBM vs. Analytic Z Velocity Profiles", computed_analytic_velocity_profiles.z, computed_lbm_velocity_profiles.z, z_velocity_error);

#ifdef VTK_EXTENSION
  Info("Writing Plot files...");
  // Plot Errors
  std::vector<VTK::CSVPlotSource> y_velocity_plot_sources;
  if(have_analytic_velocity_baseline)
    y_velocity_plot_sources.push_back(VTK::CSVPlotSource("Analytic Error", analytic_y_velocity_error, 0.0, 0.0, 0.0, 6.f, 1/*SOLID_LINE*/));
  if (have_lbm_velocity_baseline)
    y_velocity_plot_sources.push_back(VTK::CSVPlotSource("LBM VS Baseline Error", lbm_y_velocity_error, 0.0, 0.0, 1.0, 4.f, 2/*DASH_LINE*/));
  y_velocity_plot_sources.push_back(VTK::CSVPlotSource("LBM VS Analytic Error", y_velocity_error, 1.0, 0.0, 0.0, 2.f, 5/*DASH_DOT_DOT_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_full_y_velocity_error_profile_", y_velocity_plot_sources, true, *GetLogger());

  std::vector<VTK::CSVPlotSource> z_velocity_plot_sources;
  if (have_analytic_velocity_baseline)
    z_velocity_plot_sources.push_back(VTK::CSVPlotSource("Analytic Error", analytic_z_velocity_error, 0.0, 0.0, 0.0, 6.f, 1/*SOLID_LINE*/));
  if (have_lbm_velocity_baseline)
    z_velocity_plot_sources.push_back(VTK::CSVPlotSource("LBM VS Baseline Error", lbm_z_velocity_error, 0.0, 0.0, 1.0, 4.f, 2/*DASH_LINE*/));
  z_velocity_plot_sources.push_back(VTK::CSVPlotSource("LBM VS Analytic Error", z_velocity_error, 1.0, 0.0, 0.0, 2.f, 5/*DASH_DOT_DOT_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_full_z_velocity_error_profile_", z_velocity_plot_sources, true, *GetLogger());

  // Plot cross section velocities
  y_velocity_plot_sources.clear();
  y_velocity_plot_sources.push_back(VTK::CSVPlotSource("Analytic", computed_analytic_velocity_profiles.y, 0.0, 0.0, 0.0, 6.f, 1/*SOLID_LINE*/));
  if (have_lbm_velocity_baseline)
    y_velocity_plot_sources.push_back(VTK::CSVPlotSource("Baseline", baseline_lbm_velocity_profiles.y, 0.0, 0.0, 1.0, 4.f, 2/*DASH_LINE*/));
  y_velocity_plot_sources.push_back(VTK::CSVPlotSource("Computed", computed_lbm_velocity_profiles.y, 1.0, 0.0, 0.0, 2.f, 5/*DASH_DOT_DOT_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_y_velocity_profile_at_", y_velocity_plot_sources, true, *GetLogger());
  
  z_velocity_plot_sources.clear();
  z_velocity_plot_sources.push_back(VTK::CSVPlotSource("Analytic", computed_analytic_velocity_profiles.z, 0.0, 0.0, 0.0, 6.f, 1/*SOLID_LINE*/));
  if (have_lbm_velocity_baseline)
    z_velocity_plot_sources.push_back(VTK::CSVPlotSource("Baseline", baseline_lbm_velocity_profiles.z, 0.0, 0.0, 1.0, 4.f, 2/*DASH_LINE*/));
  z_velocity_plot_sources.push_back(VTK::CSVPlotSource("Computed", computed_lbm_velocity_profiles.z, 1.0, 0.0, 0.0, 2.f, 5/*DASH_DOT_DOT_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_z_velocity_profile_at_", z_velocity_plot_sources, true, *GetLogger());
#endif

  //////////////////////////////////////////////
  // Pressure, Temperature, Wall Stress Plots //
  //////////////////////////////////////////////

  //////////////////////////////////
  // Temperature Error and Timing //
  //////////////////////////////////

  LBM::CSV temperature_error;
  LBM::CSV temperature_timing;
  iterations.clear();
  for (size_t i = 0; i < m_lbm_run.out.temperature_iteration_error.size(); i++)
    iterations.push_back(i);
  temperature_error.push_back(std::pair<std::string, std::vector<float>>("iteration", iterations));
  temperature_error.push_back(std::pair<std::string, std::vector<float>>("error", m_lbm_run.out.temperature_iteration_error));
  LBM::WriteCSV(m_out_dir + m_run_name + "_temperature_iteration_error.csv", temperature_error);
  temperature_timing.push_back(std::pair<std::string, std::vector<float>>("iteration", iterations));
  temperature_timing.push_back(std::pair<std::string, std::vector<float>>("time_ms", m_lbm_run.out.temperature_iteration_timing_ms));
  LBM::WriteCSV(m_out_dir + m_run_name + "_temperature_iteration_timing.csv", temperature_timing);

#ifdef VTK_EXTENSION
  // Temperature
  std::vector<VTK::CSVPlotSource> temperature_error_source;
  temperature_error_source.push_back(VTK::CSVPlotSource("Temperature Iteration Error", temperature_error, 1.0, 0.0, 0.0, 2.f, 1/*SOLID_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_temperature_", temperature_error_source, true, *GetLogger());
  std::vector<VTK::CSVPlotSource> temperature_timing_source;
  temperature_timing_source.push_back(VTK::CSVPlotSource("Temperature Iteration Time(ms)", temperature_timing, 1.0, 0.0, 0.0, 2.f, 1/*SOLID_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_temperature_", temperature_timing_source, true, *GetLogger());
#endif

  /////////////////////////////////////////////////////////////////
  // Compare Analytic and LBM Pressure, Temperature, Wall Stress //
  /////////////////////////////////////////////////////////////////

  LBM::CSV lbm_pressure_temperature_wall_stress_baseline;
  LBM::CSV computed_lbm_pressure_temp_wall_stress_values;
  LBM::CSV computed_lbm_pressure_temp_wall_stress_error;
  bool have_lbm_pressure_temperature_wall_stress_baseline = LBM::ReadCSV(baseline_root + m_run_name + "_lbm_pressure_temperature_wall_stress_baseline.csv", lbm_pressure_temperature_wall_stress_baseline);
  // Compare and write out results
  GetPressureTemperatureWallStressValues(m_lbm_run, computed_lbm_pressure_temp_wall_stress_values, false);
  if (have_lbm_pressure_temperature_wall_stress_baseline)
  {
    m_verification_report.CompareValues(m_run_name,
                                       "Computed LBM vs. Baseline LBM Pressure, Temperature, Wall Stress",
                                        lbm_pressure_temperature_wall_stress_baseline,
                                        computed_lbm_pressure_temp_wall_stress_values,
                                        computed_lbm_pressure_temp_wall_stress_error);
  }
  else
  {
    m_verification_report.m_run_profile_comparisons[m_run_name]["MISSING Computed LBM vs. Baseline LBM Pressure, Temperature, Wall Stress Profiles"];
    status.errors.push_back("LBM pressure, temperature, wall stress baseline not found");
    status.vs_baseline_temperature = false;
  }
  // Write out csv as the new baseline, if needed
  LBM::WriteCSV(m_out_dir + m_run_name + "_lbm_pressure_temperature_wall_stress_baseline.csv", computed_lbm_pressure_temp_wall_stress_values);

  // Join the analytic and LBM CSVs into one to write
  LBM::CSV both_pressure_temp_wall_stress_values;
  for (size_t i = 0; i < computed_analytic_pressure_temperature_wall_stress_values.size(); i++)
  {
    if (i == 0)
      both_pressure_temp_wall_stress_values.push_back(computed_analytic_pressure_temperature_wall_stress_values[i]);
    else
    {
      both_pressure_temp_wall_stress_values.push_back(computed_analytic_pressure_temperature_wall_stress_values[i]);
      both_pressure_temp_wall_stress_values.back().first = "Analytic_" + computed_analytic_pressure_temperature_wall_stress_values[i].first;
      both_pressure_temp_wall_stress_values.push_back(computed_lbm_pressure_temp_wall_stress_values[i]);
      both_pressure_temp_wall_stress_values.back().first = "Computed_LBM_" + computed_analytic_pressure_temperature_wall_stress_values[i].first;
      if (have_lbm_pressure_temperature_wall_stress_baseline)
      {
        both_pressure_temp_wall_stress_values.push_back(lbm_pressure_temperature_wall_stress_baseline[i]);
        both_pressure_temp_wall_stress_values.back().first = "Baseline_LBM_" + lbm_pressure_temperature_wall_stress_baseline[i].first;
      }
    }
  }
  Info("Writing Full CSV files...");
  LBM::WriteCSV(m_out_dir + m_run_name + "_pressure_temperature_wall_stress_values.csv", both_pressure_temp_wall_stress_values);
  // Compare and write out results
  LBM::CSV computed_lbm_vs_analytic_pressure_temp_wall_stress_error;
  m_verification_report.CompareValues(m_run_name,
                                        "Computed LBM vs. Analytic Pressure, Temperature, Wall Stress Profiles",
                                        computed_analytic_pressure_temperature_wall_stress_values,
                                        computed_lbm_pressure_temp_wall_stress_values,
                                        computed_lbm_vs_analytic_pressure_temp_wall_stress_error);

#ifdef VTK_EXTENSION
  Info("Writing Plot files...");
  // Plot Errors
  std::vector<VTK::CSVPlotSource> pressure_temperature_wall_stress_plot_sources;
  if (have_analytic_pressure_temp_wall_stress_baseline)
    pressure_temperature_wall_stress_plot_sources.push_back(VTK::CSVPlotSource("Analytic Error", computed_analytic_pressure_temperature_wall_stress_error, 0.0, 0.0, 0.0, 6.f, 1/*SOLID_LINE*/));
  if (have_lbm_pressure_temperature_wall_stress_baseline)
    pressure_temperature_wall_stress_plot_sources.push_back(VTK::CSVPlotSource("LBM VS Baseline Error", computed_lbm_pressure_temp_wall_stress_error, 0.0, 0.0, 1.0, 4.f, 2/*DASH_LINE*/));
  pressure_temperature_wall_stress_plot_sources.push_back(VTK::CSVPlotSource("LBM VS Analytic Error", computed_lbm_vs_analytic_pressure_temp_wall_stress_error, 1.0, 0.0, 0.0, 2.f, 5/*DASH_DOT_DOT_LINE*/));
  //LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_pressure_temperature_wall_stress_error_profile_", pressure_temperature_wall_stress_plot_sources, true, *GetLogger());

  // Plot our values
  pressure_temperature_wall_stress_plot_sources.clear();
  pressure_temperature_wall_stress_plot_sources.push_back(VTK::CSVPlotSource("Analytic", computed_analytic_pressure_temperature_wall_stress_values, 0.0, 0.0, 0.0, 6.f, 1/*SOLID_LINE*/));
  if (have_lbm_pressure_temperature_wall_stress_baseline)
    pressure_temperature_wall_stress_plot_sources.push_back(VTK::CSVPlotSource("Baseline", lbm_pressure_temperature_wall_stress_baseline, 0.0, 0.0, 1.0, 4.f, 2/*DASH_LINE*/));
  pressure_temperature_wall_stress_plot_sources.push_back(VTK::CSVPlotSource("Computed", computed_lbm_pressure_temp_wall_stress_values, 1.0, 0.0, 0.0, 2.f, 5/*DASH_DOT_DOT_LINE*/));
  LBM::VTK::PlotCSV(m_out_dir + m_run_name + "_", pressure_temperature_wall_stress_plot_sources, false, *GetLogger());
#endif

  Info("Finished "+m_run_name);
  return;
}

void LBM::Verification::GetVelocityProfiles(LBM::Run const& lbm, LBM::VelocityProfiles& profiles, size_t num_slices)
{
  if (num_slices <= 2)
  {
    num_slices = 3;
    lbm.Info("Generating the minimum of 3 profile slices");
  }

  lbm.Info("Calculating " + to_string(num_slices) + " velocity profile slices");
  // Compute the step size down the x-axis based on the number of slices requested
  int pad = 2; // Our pad is 2 because,
  // Our verification has x0 as inactive and x1 is the inlet that does not have values
  // The same is true for the outlet and a 1 voxel pad on the other end
  size_t x_step_size = (lbm.in.dimensions[0] - 1 - pad) / (num_slices - 1);

  // Step the x and compute the profile along the y and z axis at each step

  profiles.y.resize(num_slices + 1);
  profiles.z.resize(num_slices + 1);
  auto& y_indexes = profiles.y[0];
  y_indexes.first = "y_axial_index";
  auto& z_indexes = profiles.z[0];
  z_indexes.first = "z_axial_index";

  for (size_t y = 0; y < lbm.in.dimensions[1]; y++)
    y_indexes.second.push_back((float)y);
  for (size_t z = 0; z < lbm.in.dimensions[2]; z++)
    z_indexes.second.push_back((float)z);

  size_t x_step = pad;
  for (size_t s = 0; s < num_slices; s++)
  {
    auto& y_slice = profiles.y[s + 1];
    y_slice.first = "x=" + to_string(x_step);
    auto& z_slice = profiles.z[s + 1];
    z_slice.first = "x=" + to_string(x_step);
    for (size_t y = 0; y < lbm.in.dimensions[1]; y++)
      y_slice.second.push_back(lbm.out.GetXFlow(x_step, y, (size_t)(lbm.in.dimensions[2] / 2)));
    for (size_t z = 0; z < lbm.in.dimensions[2]; z++)
      z_slice.second.push_back(lbm.out.GetXFlow(x_step, (size_t)(lbm.in.dimensions[1] / 2), z));

    x_step += x_step_size;
  }
}

void LBM::Verification::GetPressureTemperatureWallStressValues(LBM::Run const& lbm, LBM::CSV& csv, bool analytic_wall_label)
{
  csv.resize(4);
  auto& indexes = csv[0];
  indexes.first = "x_axial_index";
  auto& pressure = csv[1];
  pressure.first = "pressure";
  auto& temperature = csv[2];
  temperature.first = "temperature";
  auto& wall_stress = csv[3];
  wall_stress.first = "wall_stress";
  size_t x_center = (size_t)(lbm.in.dimensions[0] * 0.5);
  size_t y_center = (size_t)(lbm.in.dimensions[1] * 0.5);
  size_t z_center = (size_t)(lbm.in.dimensions[2] * 0.5);
  size_t y_wall = -1;
  for (size_t y = 0; y < lbm.in.dimensions[1]; y++)
  { // Walk the from y0 to ymax looking for a wall node, and get the stress from the first one we find
    int label = lbm.out.GetLabel(x_center, y, z_center);
    if ((analytic_wall_label && label==0) || (!analytic_wall_label && label>0 && label<=99 ))
    {
      y_wall = y;
      break;
    }
  }
  if (y_wall == -1)
  {
    Fatal("Unable to find a wall label for wall stress ");
    return;
  }
  for (size_t x = 0; x < lbm.in.dimensions[0]; x++)
  {
    indexes.second.push_back((float)x);
    pressure.second.push_back(lbm.out.GetPressure(x, y_center, z_center));
    temperature.second.push_back(lbm.out.GetTemperature(x, y_center, z_center));
    int l = lbm.out.GetLabel(x_center, y_wall, z_center);
    wall_stress.second.push_back(lbm.out.GetWallShearStress(x, y_wall, z_center));
  }
}

void LBM::BendVerification::RunComparison(std::string const& baseline_root, std::string const& out_dir, const LBM::Run& lbm)
{
  m_run_name = lbm.cfg.name;
  m_out_dir = out_dir;
  Logger log(m_out_dir + m_run_name + "_verification.log");
  m_lbm_run.SetLogger(log);
  m_analytic_run.SetLogger(log);
  m_lbm_run.cfg = lbm.cfg;
  m_analytic_run.cfg = lbm.cfg;
  m_lbm_run.in.grid_spacing = lbm.in.grid_spacing;
  m_analytic_run.in.grid_spacing = lbm.in.grid_spacing;
  //GenerateGeometry(m_analytic_run, radius, length);
  //GenerateGeometry(m_lbm_run, radius, length);
  Verification::RunComparison(baseline_root);
}
void LBM::BendVerification::GenerateGeometry(LBM::Run& run, float radius, float length)
{
  run.Error("Bend is currently not supported");
}
bool LBM::BendVerification::AnalyticSolution(LBM::Run& run)
{
  run.Error("Bend is currently not supported");
  // Allocate output arrays
  //run.out.Allocate();
  //run.Info("Running Bend Analytic Solution");
  return false;
}

void LBM::BrickVerification::RunComparison(std::string const& baseline_root, std::string const& out_dir, const LBM::Run& lbm)
{
  m_run_name = lbm.cfg.name;
  m_out_dir = out_dir;
  Logger log(m_out_dir + m_run_name + "_verification.log");
  m_lbm_run.SetLogger(log);
  m_analytic_run.SetLogger(log);
  m_lbm_run.cfg = lbm.cfg;
  m_analytic_run.cfg = lbm.cfg;
  m_lbm_run.in.grid_spacing = lbm.in.grid_spacing;
  m_analytic_run.in.grid_spacing = lbm.in.grid_spacing;
  GenerateGeometry(m_analytic_run, lbm.in.brick_x_length, lbm.in.brick_y_length, lbm.in.brick_z_length);
  GenerateGeometry(m_lbm_run, lbm.in.brick_x_length, lbm.in.brick_y_length, lbm.in.brick_z_length);
  Verification::RunComparison(baseline_root);
}
void LBM::BrickVerification::GenerateGeometry(LBM::Run& run, float x, float y, float z)
{
  //dimensions
  run.in.dimensions[0] = int(x / run.in.grid_spacing)+2;
  run.in.dimensions[1] = int(y / run.in.grid_spacing)+2;
  run.in.dimensions[2] = int(z / run.in.grid_spacing)+2;

  //All dimensions have to be divisible by 16 - so using modulus to pad the geometry
  /*int x_pad = run.in.dimensions[0] % 16;
  run.in.dimensions[0] += x_pad;

  int y_pad = run.in.dimensions[1] % 16;
  run.in.dimensions[1] += y_pad;

  int z_pad = run.in.dimensions[2] % 16;
  run.in.dimensions[2] += z_pad;*/

  int x_pad = 0;
  int y_pad = 0;
  int z_pad = 0;

  //labels
  for (int z = 0; z < run.in.dimensions[2]; z++)
  {
    for (int y = 0; y < run.in.dimensions[1]; y++)
    {
      for (int x = 0; x < run.in.dimensions[0]; x++)
      {
        //walls only extend to the padding, so stop at 1+pad
        if (x>0 && x<run.in.dimensions[0]-(1+x_pad) && y>0 && y<run.in.dimensions[1]-(1+y_pad) && z>0 && z<run.in.dimensions[2]-(1+z_pad))
        {
          //geometry
          if (x == 1)
          {
            //inlet
            run.in.labels.push_back(LBM::BoundaryTypes::Inlet);
          }
          else if (x == run.in.dimensions[0] - (2+x_pad)) //the outlet is at 2+pad from the end
          {
            //outlet
            run.in.labels.push_back(LBM::BoundaryTypes::Outlet);
          }
          else
          {
            //everything else is a wall
            run.in.labels.push_back(LBM::BoundaryTypes::Wall);
          }
        }
        else
        {
          //surrounding bounding box
          run.in.labels.push_back(LBM::BoundaryTypes::Inactive);
        }
      }
    }
  }
}

bool LBM::BrickVerification::AnalyticSolution(LBM::Run& run)
{
  if(run.cfg.inlet_boundary_condition != BoundaryCondition::Pressure ||
     run.cfg.outlet_boundary_condition != BoundaryCondition::Pressure)
  {
    run.Error("Currently the brick analytic solution does not support flow boundary conditions");
    return false;
  }
  // Allocate output arrays
  run.out.Allocate();
  run.Info("Running Brick Analytic Solution");
  float spacing = run.in.grid_spacing;
  size_t nx = run.in.dimensions[0];
  size_t ny = run.in.dimensions[1];
  size_t nz = run.in.dimensions[2];
  float viscosity = run.cfg.viscosity;
  float density = run.cfg.density;
  float Pin = run.cfg.inlet_pressure;
  float Pout = run.cfg.outlet_pressure;
  size_t num_cells = run.out.num_cells;

  float delP = Pin - Pout;
  float length = 0.08;//(nx - 2) * spacing;
  float length_y = 0.02;//(ny)*spacing;
  float length_z = 0.01;//(nz)*spacing;
  float friction_coeff = delP / (2 * viscosity * density * length);

  size_t iCell = 0;
  for (size_t z = 0; z < nz; z++)
  {
    for (size_t y = 0; y < ny; y++)
    {
      for (size_t x = 0; x < nx; x++)
      {
        int label = (int)run.in.labels[iCell];
        run.in.source_labels.push_back(label);
        run.out.labels[iCell] = label;
        if (label >= 0 && label <= 300)
        {
          size_t N = 15;
          float series_sum_u = 0.f;
          float series_sum_q = 0.f;

          for (size_t i = 1; i <= N; i++)
          {
            float beta = M_PI * (2 * i - 1) / (length_y);
            float factor_u = 1 / std::pow((2 * i - 1), 3);
            float factor_q = 1 / std::pow((2 * i - 1), 5);
            float tmp_u = (sinh(beta * z * spacing) + sinh(beta * (length_z - z * spacing))) * sin(beta * y * spacing) / sinh(beta * length_z);
            float tmp_q = (cosh(beta * length_z) - 1) / sinh(beta * length_z);
            series_sum_u += factor_u * tmp_u;
            series_sum_q += factor_q * tmp_q;
          }

          run.out.pressure_flows_and_stress[iCell] = Pout * x / nx; // linear pressure drop
          run.out.pressure_flows_and_stress[num_cells + iCell] = friction_coeff * (y * spacing * (length_y - y * spacing) - 8 * std::pow(length_y, 2) * series_sum_u / std::pow(M_PI, 3)); //axial velocity
          run.out.pressure_flows_and_stress[(num_cells * 2) + iCell] = friction_coeff * (std::pow(length_y, 3) * length_z / 6 - 32 * std::pow(length_y, 4) * series_sum_q / std::pow(M_PI, 5)); // volume flow-rate
          run.out.pressure_flows_and_stress[(num_cells * 3) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 4) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 5) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 6) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 7) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 8) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 9) + iCell] = 0.f;  // place holder for wall shear stress
          run.out.wall_shear_stress[iCell] = 0.f; // wall shear stress
          run.out.wall_shear_stress[num_cells + iCell] = 0.f;
        }
        else
        {
          run.out.pressure_flows_and_stress[iCell] = 0.f;
          run.out.pressure_flows_and_stress[num_cells + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 2) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 3) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 4) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 5) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 6) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 7) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 8) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 9) + iCell] = 0.f;
          run.out.wall_shear_stress[iCell] = 0.f;
          run.out.wall_shear_stress[num_cells + iCell] = 0.f;
        }
        iCell++;
      }
    }
  }
  return true;
}

void LBM::CylinderVerification::RunComparison(std::string const& baseline_root, std::string const& out_dir, const LBM::Run& lbm)
{
  m_run_name = lbm.cfg.name;
  m_out_dir = out_dir;
  Logger log(m_out_dir + m_run_name + "_verification.log");
  m_lbm_run.SetLogger(log);
  m_analytic_run.SetLogger(log);
  m_lbm_run.cfg = lbm.cfg;
  m_analytic_run.cfg = lbm.cfg;
  m_lbm_run.in.grid_spacing = lbm.in.grid_spacing;
  m_analytic_run.in.grid_spacing = lbm.in.grid_spacing;
  GenerateGeometry(m_analytic_run, lbm.in.cylinder_radius, lbm.in.cylinder_length);
  GenerateGeometry(m_lbm_run, lbm.in.cylinder_radius, lbm.in.cylinder_length);
  Verification::RunComparison(baseline_root);
}
void LBM::CylinderVerification::GenerateGeometry(LBM::Run& run, float radius, float length)
{
  //dimensions
  run.in.dimensions[0] = int(length / run.in.grid_spacing) + 2;
  run.in.dimensions[1] = int(radius*2 / run.in.grid_spacing) + 2;
  run.in.dimensions[2] = int(radius*2 / run.in.grid_spacing) + 2;

  //get the radius
  int radius_vox = (int)(run.in.dimensions[1] - 2) * 0.5f;

  //All dimensions have to be divisible by 16 - so using modulus to pad the geometry
  /*int x_pad = run.in.dimensions[0] % 16;
  run.in.dimensions[0] += x_pad;

  int y_pad = run.in.dimensions[1] % 16;
  run.in.dimensions[1] += y_pad;

  int z_pad = run.in.dimensions[2] % 16;
  run.in.dimensions[2] += z_pad;*/

  int x_pad = 0;
  int y_pad = 0;
  int z_pad = 0;

  int y_center = radius_vox + 1;
  int z_center = radius_vox + 1;

  //labels
  for (int z = 0; z < run.in.dimensions[2]; z++)
  {
    for (int y = 0; y < run.in.dimensions[1]; y++)
    {
      for (int x = 0; x < run.in.dimensions[0]; x++)
      {
        //walls only extend to the padding, so stop at 1+pad
        int distance = std::pow(std::pow(y-y_center, 2) + std::pow(z-z_center, 2), 0.5);
        if (x>0 && x<run.in.dimensions[0] - (1 + x_pad) && distance<radius_vox)
        {
          //geometry
          if (x == 1)
          {
            //inlet
            run.in.labels.push_back(LBM::BoundaryTypes::Inlet);
          }
          else if (x == run.in.dimensions[0] - (2 + x_pad)) //the outlet is at 2+pad from the end
          {
            //outlet
            run.in.labels.push_back(LBM::BoundaryTypes::Outlet);
          }
          else
          {
            //everything else is a wall
            run.in.labels.push_back(LBM::BoundaryTypes::Wall);
          }
        }
        else
        {
          //surrounding bounding box
          run.in.labels.push_back(LBM::BoundaryTypes::Inactive);
        }
      }
    }
  }

}
bool LBM::CylinderVerification::AnalyticSolution(LBM::Run& run)
{
  if(run.cfg.inlet_boundary_condition != BoundaryCondition::Pressure ||
     run.cfg.outlet_boundary_condition != BoundaryCondition::Pressure)
  {
    run.Error("Currently the cylinder analytic solution does not support flow boundary conditions");
    return false;
  }
  // Allocate output arrays
  run.out.Allocate();
  run.Info("Running Cylinder Analytic Solution");
  // Analytical solution of steady-state Navier-Stokes equation in an uniform cylinder of finite length
  // Axial velocity profile (Poiseuille's profile), volume flow rate, and wall shear stress
  float spacing = run.in.grid_spacing;
  size_t nx = run.in.dimensions[0];
  size_t ny = run.in.dimensions[1];
  size_t nz = run.in.dimensions[2];
  float viscosity = run.cfg.viscosity;
  float density = run.cfg.density;
  float Pin = run.cfg.inlet_pressure;
  float Pout = run.cfg.outlet_pressure;
  size_t num_cells = run.out.num_cells;

  float delP = Pin - Pout;
  float length = 0.1;//(nx - 2) * spacing;
  float rad =  0.01;//ny / 2 * spacing;
  float rad2 = std::pow(rad, 2);
  float rad3 = std::pow(rad, 3);
  float rad4 = std::pow(rad, 4);
  float Geometric_Resistance = delP * rad2 / (4 * viscosity * density * length);

  size_t iCell = 0;
  for (size_t z = 0; z < nz; z++)
  {
    for (size_t y = 0; y < ny; y++)
    {
      for (size_t x = 0; x < nx; x++)
      {
        int label = (int)run.in.labels[iCell];
        run.in.source_labels.push_back(label);
        run.out.labels[iCell] = label;
        if (label >= 0 && label <= 300)
        {
          run.out.pressure_flows_and_stress[iCell] = Pout * x / nx; // linear pressure drop
          run.out.pressure_flows_and_stress[num_cells + iCell] = Geometric_Resistance * (1 - (std::pow((((z + 0.5) * spacing - rad) / rad), 2) + std::pow((((y + 0.5) * spacing - rad) / rad), 2))); //axial velocity
          run.out.pressure_flows_and_stress[(num_cells * 2) + iCell] = M_PI * rad4 * delP / (8 * viscosity * density * length); // volume flow rate
          run.out.pressure_flows_and_stress[(num_cells * 3) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 4) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 5) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 6) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 7) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 8) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 9) + iCell] = 0.f;
          run.out.wall_shear_stress[iCell] = 4 * viscosity * density * run.out.pressure_flows_and_stress[(num_cells * 2) + iCell] / (M_PI * rad3); // wall shear stress
          run.out.wall_shear_stress[num_cells + iCell] = 0.f;
        }
        else
        {
          run.out.pressure_flows_and_stress[iCell] = 0.f;
          run.out.pressure_flows_and_stress[num_cells + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 2) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 3) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 4) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 5) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 6) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 7) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 8) + iCell] = 0.f;
          run.out.pressure_flows_and_stress[(num_cells * 9) + iCell] = 0.f;
          run.out.wall_shear_stress[iCell] = 0.f;
          run.out.wall_shear_stress[num_cells + iCell] = 0.f;
        }
        iCell++;
      }
    }
  }
  return true;
}

bool LBM::VerificationReport::CompareProfiles(ComparisonTypes type, std::string const& run_name, std::string const& table_name, const LBM::CSV& baseline, const LBM::CSV& computed, LBM::CSV& error)
{
  bool pass = true;

  // Create a CSV for the error for plotting
  float x_value;
  std::vector<float> x_index;
  std::vector<float> average_error;
  std::vector<float> average_error_fraction;

  auto& result = m_run_profile_comparisons[run_name][table_name];
  auto& status = m_run_status[run_name];
  if (baseline[0].first != computed[0].first)
  {
    status.errors.push_back("Computed CSV abscissa does not match basline abscissa");
    return false;
  }

  float max_error = -std::numeric_limits<float>::max();
  switch (type)
  {
  case ComparisonTypes::AnalyticVelocity:
    max_error = status.max_analytic_velocity_error;
    break;
  case ComparisonTypes::VS_BaselineVelocity:
    max_error = status.max_vs_baseline_velocity_error;
    break;
  case ComparisonTypes::VS_AnalyticVelocity:
    max_error = status.max_vs_analytic_velocity_error;
    break;
  }

  for (size_t b=1; b<baseline.size(); b++)
  {
    LBM::ProfileComparison pc;

    // Parse out the x value from the string name
    pc.name = baseline[b].first;
    x_value = std::stof(pc.name.substr(pc.name.find("x=") + 2));
    x_index.push_back(x_value);
    const std::vector<float>* b_vector = &baseline[b].second;
    const std::vector<float>* c_vector = nullptr;
    for (size_t c=1; c< computed.size(); c++)
    {
      if (computed[c].first == pc.name)
      {
        c_vector = &computed[c].second;
        break;
      }
    }
    if (c_vector == nullptr)
    {
      pass = false;
      status.errors.push_back("Unable to find "+ pc.name +" in computed results");
      continue;
    }
    if (c_vector->size() != b_vector->size())
    {
      pass = false;
      status.errors.push_back("Vectors for header "+ pc.name +" are not the same size");
      continue;
    }

    float avg_vel = 0;
    pc.number_voxels = b_vector->size();
    pc.number_voxel_failures = 0;
    pc.average_error = 0;
    for (size_t i=0; i<b_vector->size(); i++)
    {
      float b_vel = (*b_vector)[i];
      float c_vel = (*c_vector)[i];
      double diff = LBM::PercentDifference(b_vel, c_vel);
      if (diff > percrent_difference)
      {
        pass = false;
        pc.number_voxel_failures++;
      }
      //pc.least_mean_squared_error += std::pow(b_vel - c_vel, 2);
      pc.average_error += (b_vel - c_vel);
      avg_vel += b_vel;
    }
    avg_vel /= pc.number_voxels;
    //pc.least_mean_squared_error = sqrt(pc.least_mean_squared_error);
    pc.average_error /= pc.number_voxels;
    average_error.push_back(pc.average_error);
    pc.average_error_fraction = pc.average_error / avg_vel;
    average_error_fraction.push_back(pc.average_error_fraction);
    if (pc.average_error_fraction > max_error)
      max_error = pc.average_error_fraction;

    result.push_back(pc);
  }

  // Construct our CSV from our error
  error.clear();
  error.push_back(std::pair <std::string, std::vector<float>>("x_index", x_index));
  error.push_back(std::pair <std::string, std::vector<float>>("average_error", average_error));
  error.push_back(std::pair <std::string, std::vector<float>>("average_error_fraction", average_error_fraction));

  // Add to our run status
  switch (type)
  {
  case ComparisonTypes::AnalyticVelocity:
    status.analytic_velocity = pass;
    status.max_analytic_velocity_error = max_error;
    break;
  case ComparisonTypes::VS_BaselineVelocity:
    status.vs_baseline_velocity = pass;
    status.max_vs_baseline_velocity_error = max_error;
    break;
  case ComparisonTypes::VS_AnalyticVelocity:
    //status.vs_analytic_velocity = pass;
    status.max_vs_analytic_velocity_error = max_error;
    // Judging LBM vs Analytic to be based only the middle slice
    size_t mid = average_error_fraction.size() * 0.75;
    float mid_lbm_error = average_error_fraction[mid];
    if (mid_lbm_error > status.mid_vs_analytic_velocity_error)
      status.mid_vs_analytic_velocity_error = mid_lbm_error;
    if(status.mid_vs_analytic_velocity_error*100 > 10)
      status.vs_analytic_velocity = false;
    break;
  }

  return pass;
}

bool LBM::VerificationReport::CompareValues(std::string const& run_name, std::string const& table_name, const LBM::CSV& baseline, const LBM::CSV& computed, LBM::CSV& error)
{
  return false;
}

bool LBM::VerificationReport::RecordRuntime(std::string const& run_name, double runtime_s)
{
  if (m_run_status.find(run_name) == m_run_status.end())
    return false;
  auto& status = m_run_status[run_name];
  status.total_time_to_run_s = runtime_s;
  return true;
}

void LBM::VerificationReport::WriteReport(std::string const& to_dir)
{
  std::string green = "bgcolor=\"#00FF00\"";
  std::string yellow = "bgcolor=\"#FFFF00\"";
  std::string red = "bgcolor=\"#FF0000\"";

  std::string* color;
  std::string  status;
  float totalTime_s=0;

  std::ofstream rpt;
  rpt.open(to_dir + "LBM_Verification_Report.html", std::ofstream::out | std::ofstream::trunc);
  if (!rpt.good())
  {
    Error("Unable to write to directory " + to_dir);
    return;
  }

  rpt << "<html>\n";
  rpt << "<head><title>LBM Verification Report</title></head>\n";
  rpt << "<body><br>\n";
  rpt << "<center>\n";

  rpt << "<table border=\"1\">\n";
  rpt << "<tr>\n";
  rpt << "<th>  Run  </th>\n";
  rpt << "<th>  Total Test time </th>\n";
  rpt << "<th>  LBM Compute time </th>\n";
  rpt << "<th>  Num LBM Iterations </th>\n";
  rpt << "<th>  Num Nodes </th>\n";
  rpt << "<th>  MLUPS  </th>\n";
  rpt << "<th>  Velocity Baseline Analytic vs. Computed Analytic  </th>\n";
  rpt << "<th>  Max Avg Velocity Error Fraction Baseline LBM vs. Computed LBM  </th>\n";
  rpt << "<th>  Average Velocity Error Fraction Computed LBM vs Computed Analytic  </th>\n";
  rpt << "<th>  Execution Errors  </th>\n";
  rpt << "</tr>\n\n";

  for (auto& run : m_run_status)
  {
    rpt << "<tr>\n";
    rpt << "<td>" + run.first + "</td>\n";
    rpt << "<td><center>" + to_string(run.second.total_time_to_run_s) + "s </center></td>\n";
    rpt << "<td><center>" + to_string(run.second.lbm_compute_time_s) + "s </center></td>\n";
    rpt << "<td><center>" + (run.second.num_fluid_iterations == 0 ? "N/A" : to_string(run.second.num_fluid_iterations)) + " </center></td>\n";
    rpt << "<td><center>" + to_string(run.second.num_nodes) + " </center></td>\n";
    rpt << "<td><center>" + (run.second.num_fluid_iterations == 0 ? "N/A" : to_string( (run.second.num_nodes / 1000000.0) / (run.second.lbm_compute_time_s / run.second.num_fluid_iterations))) + " </center></td>\n";

    color = run.second.analytic_velocity ? &green : &red;
    status = run.second.analytic_velocity ? " (Pass)" : " (Fail)";
    rpt << "<td " + (*color) + " ><center> "+to_string(run.second.max_analytic_velocity_error)+status+" </center></td>\n";
    color = run.second.vs_baseline_velocity ? &green : &red;
    status = run.second.vs_baseline_velocity ? " (Pass)" : " (Fail)";
    rpt << "<td " + (*color) + " ><center> " + to_string(run.second.max_vs_baseline_velocity_error)+status+" </center></td>\n";
    color = run.second.vs_analytic_velocity ? &green : &red;
    status = run.second.vs_analytic_velocity ? " (Pass)" : " (Fail)";
    rpt << "<td " + (*color) + " ><center> Mid : " + to_string(run.second.mid_vs_analytic_velocity_error)+"<br>" 
      + "Max : " + to_string(run.second.max_vs_analytic_velocity_error)+status+" </center></td>\n";
    if(run.second.errors.empty())
      rpt << "<td><center>None</center></td>\n";
    else
    {
      rpt << "<td><center>\n";
      for (std::string error : run.second.errors)
      {
        rpt << error << "<br>\n";
      }
      rpt << "</center></td>\n";
    }
    totalTime_s += run.second.total_time_to_run_s;
    rpt << "</tr>\n";
  }
  rpt << "</table>\n";
  rpt << "</center>\n";
  rpt << "<br><br><br>\n";

  int mins = (int)totalTime_s / 60;
  int secs = (int)totalTime_s % 60;
  rpt << "<b> Total Execution Time : </b>" + std::to_string(mins) + ":" + std::to_string(secs) + " (min:s) <br> \n";


  rpt << "</body>\n";
  rpt << "</html>\n";
  rpt.close();

  // Write Comparison Tables to run output folders
  for (auto& run : m_run_profile_comparisons)
  {
    std::ofstream run_rpt;
    run_rpt.open(to_dir + "/" + run.first + "/" + run.first + "_comparison_report.html", std::ofstream::out | std::ofstream::trunc);
    if (!run_rpt.good())
    {
      Error("Unable to write to directory " + to_dir + "/" + run.first);
      continue;
    }

    run_rpt << "<html>\n";
    run_rpt << "<head><title>"+run.first+" Verification Report</title></head>\n";
    run_rpt << "<body><br>\n";
    run_rpt << "<center>\n";


    for (auto& table : run.second)
    {
      run_rpt << "<table border=\"1\">\n";
      run_rpt << "<tr>\n";
      run_rpt << "<th>  " + table.first + "  </th>\n";
      run_rpt << "<th>  # of Voxels  </th>\n";
      run_rpt << "<th>  # Voxel Errors  </th>\n";
      run_rpt << "<th>  Avg. Velocity Error Fraction </th>\n";
      run_rpt << "</tr>\n\n";

      for (auto& vpc : table.second)
      {
        color = &green;
        if (vpc.number_voxel_failures > 0)
          color = &red;
        else if (vpc.average_error_fraction > 0.001)
          color = &yellow;

        run_rpt << "<tr " + (*color) + " >\n";
        run_rpt << "<td><center>  " + vpc.name + "  </center></td>\n";

        run_rpt << "<td><center>" + to_string(vpc.number_voxels) + "</center></td>\n";
        run_rpt << "<td><center>" + to_string(vpc.number_voxel_failures) + "</center></td>\n";
        run_rpt << "<td><center>" + to_string(vpc.average_error_fraction) + "</center></td>\n";
        run_rpt << "</tr>\n\n";
      }
      run_rpt << "</table>\n";
      run_rpt << "<br><br>\n";
    }

    run_rpt << "</center>\n";
    run_rpt << "</body>\n";
    run_rpt << "</html>\n";
    run_rpt.close();
  }
}
