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

#include "LBM.h"
#include "FileUtils.h"
#include <iostream>
#include <fstream>
#include <sstream>
#define _USE_MATH_DEFINES
#include <cmath>
#include <math.h>
#include <cstring>
#include<Eigen/Dense>

#include "utils/TimingProfile.h"

#define SAFE_DELETE_ARRAY(ary) \
  if(ary!=nullptr) \
  { \
    delete[] ary; \
    ary = nullptr; \
  }

   // CUDA Methods
extern "C" int microLBM(LBM::Run& run);

LBM::Input::Input()
{
  origin[0] = 0;
  origin[1] = 0;
  origin[2] = 0;
  dimensions[0] = 0;
  dimensions[1] = 0;
  dimensions[2] = 0;
  grid_spacing = 1;
}

LBM::Output::Output(Run const& r) : run(r)
{

}
LBM::Output::~Output()
{
  Deallocate();
}
bool LBM::Output::Allocate()
{
  Deallocate();

  num_cells = run.in.dimensions[0] * run.in.dimensions[1] * run.in.dimensions[2];
  if (num_cells <= 0)
  {
    Error("Input has no dimensions");
    return false;
  }
  num_wall_cells = 0;
  num_inflow_cells = 0;
  num_outflow_cells = 0;
  num_inactive_cells = 0;
  num_interior_cells = 0;
  num_unknown_cells = 0;

  Q_num_cell_values = num_cells * (nQ + nS);
  T_num_cell_values = num_cells * nQ;
  num_boundary_values = nQ * nBound;
  nPDF = num_cells * 19;

  inlet_velocity = 0;
  inlet_pressure = 0;
  outlet_pressure = 0;
  inlet_area = 0;
  outlet_velocity = 0;
  outlet_area = 0;
  average_diameter = 0;
  axial_length = 0;
  average_area = 0;
  mass_flow = 0;
  pressure_drop = 0;
  volumetric_flow_rate = 0;


  labels = new int[num_cells];
  pressure_flows_and_stress = new float[Q_num_cell_values];
  wall_shear_stress = new float[num_cells * 2];
  temperature = new float[T_num_cell_values];
  fPDF = new float[nPDF];
  fnPDF = new float[nPDF];
  boundary_conditions = new float[num_boundary_values];
  previous_results = new float[num_cells * 3];
  memset(boundary_conditions, 0, num_boundary_values * sizeof(float));

  return true;
}
void LBM::Output::Deallocate()
{
  SAFE_DELETE_ARRAY(labels);
  SAFE_DELETE_ARRAY(pressure_flows_and_stress);
  SAFE_DELETE_ARRAY(wall_shear_stress);
  SAFE_DELETE_ARRAY(temperature);
  SAFE_DELETE_ARRAY(fPDF);
  SAFE_DELETE_ARRAY(fnPDF);
  SAFE_DELETE_ARRAY(boundary_conditions);
  SAFE_DELETE_ARRAY(previous_results);
}

bool LBM::Output::CountCellTypes()
{
  if (num_cells <= 0 || labels == nullptr)
    return false;

  // Count up label types
  num_wall_cells = 0;
  num_inflow_cells = 0;
  num_outflow_cells = 0;
  num_inactive_cells = 0;
  num_interior_cells = 0;
  num_unknown_cells = 0;
  for (size_t i = 0; i < num_cells; ++i)
  {
    int label = labels[i];

    if (label == -1) // inactive - outside of the region of interest
    {
      num_inactive_cells++;
    }
    else if (label == 0) // interior of the region of interest
    {
      num_interior_cells++;
    }
    else if (label > 0 && label <= 99) // wall boundary condition area
    {
      num_wall_cells++;
    }
    else if (label > 100 && label <= 199) // inlet boundary condition area
    {
      num_inflow_cells++;
      /*if(m_gridXValues[iNodes] < xInflowStart)
          xInflowStart = gridXValues[iNodes];*/
    }
    else if (label > 200 && label <= 299) // outlet boundary condition area
    {
      num_outflow_cells++;
      /*if(m_gridXValues[iNodes] > xOutflowEnd)
          xOutflowEnd = gridXValues[iNodes];*/
    }
    else
    {
      num_unknown_cells++;
    }
  }
  Info("Number of Nodes: " + to_string(num_cells));
  Info("Number of Inactive Nodes: " + to_string(num_inactive_cells));
  Info("Number of Interior Nodes: " + to_string(num_interior_cells));
  Info("Number of Wall Nodes: " + to_string(num_wall_cells));
  Info("Number of Inflow Nodes: " + to_string(num_inflow_cells));
  Info("Number of Outflow Nodes: " + to_string(num_outflow_cells));
  Info("Number of Unknown Nodes: " + to_string(num_unknown_cells));

  return true;
}

size_t LBM::Output::Index(size_t x, size_t y, size_t z) const
{
  return (run.in.dimensions[1] * run.in.dimensions[0] * z) + (run.in.dimensions[0] * y) + x;
}

LBM::Run::Run(Logger* logger) : Loggable(logger), out(*this)
{

}
LBM::Run::Run(std::string const& logfile) : Loggable(logfile), out(*this)
{

}
LBM::Run::~Run()
{
  out.Deallocate();
}
bool LBM::Run::Compute()
{
  out.compute_time_s = 0;
  TimingProfile timer;

  if (in.labels.empty() && in.source_labels.empty())
  {
    Error("Input has no labels");
    return false;
  }
  if (in.labels.empty() && !in.source_labels.empty() && in.source_to_lbm_label_map.empty())
  {
    Error("When providing source labels, you must provide a source to LBM label mapping");
    return false;
  }
  // Check that the number of labels matches the dimensions
  size_t num_labels = in.labels.empty() ? in.source_labels.size() : in.labels.size();
  size_t num_cells = in.dimensions[0] * in.dimensions[1] * in.dimensions[2];
  if (num_labels != num_cells)
  {
    Error("Number of labels does not match the provided dimensions");
    return false;
  }

  if (in.grid_spacing <= 0)
  {
    Error("Input has no grid spacing");
    return false;
  }

  // Allocate our output
  if (!out.Allocate())
  {
    return false;
  }

  std::stringstream ss;
  Info("Processing inputs...");
  Info(ss << cfg);

  // Process labels
  size_t num_inactive = 0;
  size_t num_wall = 0;
  size_t num_inlet = 0;
  size_t num_outlet = 0;
  // Pull the initial labels, CUDA can change these
  if (!in.labels.empty())
  {
    Info("Processing provided LBM labels");
    for (size_t i = 0; i < out.num_cells; i++)
    {
      BoundaryTypes t = in.labels[i];
      out.labels[i] = (int)t;
      // Count what we have
      switch (t)
      {
      case BoundaryTypes::Inactive:
        num_inactive++; break;
      case BoundaryTypes::Wall:
        num_wall++; break;
      case BoundaryTypes::Inlet:
        num_inlet++; break;
      case BoundaryTypes::Outlet:
        num_outlet++; break;
      }
    }
  }
  else
  {
    Info("Translating source labels to LBM labels");
    for (size_t i = 0; i < out.num_cells; i++)
    {
      int source_label = in.source_labels[i];
      auto value = in.source_to_lbm_label_map.find(source_label);
      if (value == in.source_to_lbm_label_map.end())
      {
        Error("A mapping was not provided for source label " + source_label);
        return false;
      }
      out.labels[i] = (int)value->second;
      // Count what we have
      switch (value->second)
      {
      case BoundaryTypes::Inactive:
        num_inactive++; break;
      case BoundaryTypes::Wall:
        num_wall++; break;
      case BoundaryTypes::Inlet:
        num_inlet++; break;
      case BoundaryTypes::Outlet:
        num_outlet++; break;
      }
    }
  }
  Info("Labeled mesh with " + to_string(num_labels) + " labels");
  Info("  Inactive cells  : " + to_string(num_inactive));
  Info("  Wall Cells      : " + to_string(num_wall));
  Info("  Inlet Cells     : " + to_string(num_inlet));
  Info("  Outlet Cells    : " + to_string(num_outlet));

  int status = GetNormals();
  if (status)
  {
    Info("Successfully calculated the normals to the boundary.");
  }
  else
  {
    Error("Error occurred while calculating the normals.");
    return false;
  }

  Info("Find Normals Boundary Codes:");
  Info("Grid Size (x, y, z) = (" + to_string(in.dimensions[0]) + ","
    + to_string(in.dimensions[1]) + ","
    + to_string(in.dimensions[2]) + ").");
  out.CountCellTypes();

  // Calculate average area, diameter, and axial length
  out.inlet_area = out.num_inflow_cells * std::pow(in.grid_spacing, 2);
  out.outlet_area = out.num_outflow_cells * std::pow(in.grid_spacing, 2);
  out.average_area = 0.5f * (out.inlet_area + out.outlet_area);
  out.average_diameter = 2.0f * std::pow(out.average_area / float(M_PI), 0.5f);
  out.axial_length = (in.dimensions[0] - 2) * in.grid_spacing;
  out.inlet_pressure = cfg.atmospheric_pressure;

  // Organize boundary condition data specified for the continuum mechanics components for the default case
  // Given pressure drop, calculate mass and volumetric flow-rates
  if (cfg.inlet_boundary_condition == BoundaryCondition::Pressure && cfg.outlet_boundary_condition == BoundaryCondition::Pressure)
  {
    if (cfg.inlet_pressure <= cfg.outlet_pressure)
    {
      Error("!!! Inlet pressure is lower than outlet pressure !!!");
    }
    out.pressure_drop = cfg.inlet_pressure - cfg.outlet_pressure;
    out.mass_flow = std::pow(out.average_area, 2) * out.pressure_drop / (8.f * float(M_PI) * cfg.viscosity) / out.axial_length;
    out.volumetric_flow_rate = out.mass_flow / cfg.density;
    // Calculate outlet pressure and velocities
    out.outlet_pressure = out.inlet_pressure - out.pressure_drop;
//    out.inlet_velocity = out.volumetric_flow_rate / out.inlet_area;
//    out.outlet_velocity = out.volumetric_flow_rate / out.outlet_area;
    out.inlet_velocity = cfg.target_flow_rate * 1.666667e-05f / out.inlet_area;//.volumetric_flow_rate / out.inlet_area;
    out.outlet_velocity = cfg.target_flow_rate * 1.666667e-05f / out.outlet_area;//.volumetric_flow_rate / out.outlet_area;
  }
  // Given constant volumetric flow rate, calculate mass flow rate and pressure drop.
  else if (cfg.inlet_boundary_condition == BoundaryCondition::Flow && cfg.outlet_boundary_condition == BoundaryCondition::Flow)
  {
    if (cfg.inlet_volumetric_flow_rate != cfg.outlet_volumetric_flow_rate)
    {
      Error("!!! Inlet and outlet flows are NOT the same !!!");
    }
    out.volumetric_flow_rate = cfg.outlet_volumetric_flow_rate;
    out.mass_flow = out.volumetric_flow_rate * cfg.density;
    out.pressure_drop = out.mass_flow * 8.f * float(M_PI) * cfg.viscosity * out.axial_length / std::pow(out.average_area, 2);
    // Calculate outlet pressure and velocities
    out.outlet_pressure = out.inlet_pressure - out.pressure_drop;
    out.inlet_velocity = out.volumetric_flow_rate / out.inlet_area;
    out.outlet_velocity = out.volumetric_flow_rate / out.outlet_area;
  }
  else
  {
    Error("Unsupported Boundary Condition Combination.");
    return false;
  }

  float lattice_speed = 0.577350269f;  // LBM model sound speed (in non-dimensional lattice units) D3Q19:1/sqrt(3)

  // Constant conversion of factors (reference quantities) for non-dimensionalization
  // This conversion algorithm/strategy is coded for density, speed of sound & grid spacing (via the cfg file)
  // We will compute other cfg values based on those numbers
  // In the future, we will allow users to enter other cfg parametes as the base, and derive values from those
  // ex. (pressure, speed of sound, and grid spacing) or (viscosity, speed of sound and grid spacing)

  cfg.temporal_step = cfg.unitless_viscosity * std::pow(in.grid_spacing,2) /cfg.viscosity;// / (in.grid_spacing * lattice_speed / cfg.speed_of_sound);
  cfg.relaxation_time = cfg.viscosity / std::pow(cfg.speed_of_sound, 2) + 0.5f;
  cfg.relaxation_time_temperature = cfg.thermal_coefficient / std::pow(cfg.speed_of_sound, 2) + 0.5f;
  cfg.reynolds_number = out.average_diameter * out.outlet_velocity / cfg.viscosity;
  cfg.mach_number = out.outlet_velocity / cfg.speed_of_sound;

  // Conversion factors for converting physical quantities to dimensionless simulation quantities and vice-versa
  float conv_length = 0;
  float conv_time = 0;
  float conv_density = 0;
  float conv_viscosity = 0;
  float conv_velocity = 0;
  float conv_pressure = 0;

  conv_density = cfg.density;
  conv_length = in.grid_spacing;
  conv_time = cfg.temporal_step;
  conv_velocity = conv_length / conv_time; // Change this to change the time-step
  conv_pressure = conv_density * std::pow(conv_velocity, 2);
  conv_viscosity = std::pow(conv_length, 2) / conv_time;

  in.density = cfg.density / conv_density;
  in.spatial_step = in.grid_spacing / conv_length;
  in.temporal_step = cfg.temporal_step / conv_time;
  in.speed_of_sound = lattice_speed;
  in.system_length = out.average_diameter / conv_length;
  in.viscosity = cfg.viscosity / conv_viscosity;
  in.thermal_coefficient = cfg.thermal_coefficient;

  in.system_velocity = out.outlet_velocity / conv_velocity;
  in.reynolds_number = in.system_velocity * in.system_length / in.viscosity;
  in.reynolds_grid_number = in.system_velocity * in.spatial_step / in.viscosity;
  in.relaxation_time = in.viscosity / std::pow(lattice_speed, 2) + 0.5f;
  in.relaxation_time_temperature = in.thermal_coefficient / std::pow(lattice_speed, 2) + 0.5f;
  in.mach_number = in.system_velocity / in.speed_of_sound;
  in.knudsen_number = in.mach_number / in.reynolds_number;

  Info("=================== Computational Parameters ====================");
  Info("Computational Parameters:");
  Info("Grid spacing                 =    " + to_string(in.grid_spacing));
  Info("Time discretization          =    " + to_string(cfg.temporal_step));
  Info("Typical System Velocity      =    " + to_string(out.outlet_velocity));
  Info("Kinematic Viscosity          =    " + to_string(cfg.viscosity));
  Info("Dimensionless Viscosity      =    " + to_string(cfg.unitless_viscosity));
  Info("Macroscopic Density          =    " + to_string(cfg.density));
  Info("BGK Relaxation Time          =    " + to_string(cfg.relaxation_time));
  Info("BGK Lattice Relaxation Time  =    " + to_string(in.relaxation_time));
  Info("Typical System Length        =    " + to_string(out.average_diameter));

  Info("=================== Dimensionless Numbers =======================");
  Info("Reynolds Number       =    " + to_string(in.reynolds_number));
  Info("Grid Reynolds Number  =    " + to_string(in.reynolds_grid_number));
  Info("Mach Number           =    " + to_string(in.mach_number));
  Info("Knudsen Number        =    " + to_string(in.knudsen_number));

  Info("=================== Conversion Factors ==========================");
  Info("Characteristic Length         =    " + to_string(conv_length));
  Info("Characteristic Density        =    " + to_string(conv_density));
  Info("Characteristic Time           =    " + to_string(conv_time));
  Info("Characteristic Velocity       =    " + to_string(conv_velocity));
  Info("Characteristic Viscosity      =    " + to_string(conv_viscosity));
  Info("Characteristic Pressure       =    " + to_string(conv_pressure));

  Info("================ Macroscopic Boundary Conditions ================");
  // Imposed pressure drop
  if (cfg.inlet_boundary_condition == BoundaryCondition::Pressure && cfg.outlet_boundary_condition == BoundaryCondition::Pressure)
  {
    Info("Imposed pressure drop         =    " + to_string(out.pressure_drop));
    Info("Initial mass flow rate        =    " + to_string(out.mass_flow));
    Info("Initial volumetric flow rate  =    " + to_string(out.volumetric_flow_rate));
  }
  // Constant volumetric flow rate
  else if (cfg.inlet_boundary_condition == BoundaryCondition::Flow && cfg.outlet_boundary_condition == BoundaryCondition::Flow)
  {
    Info("Initial pressure drop         =    " + to_string(out.pressure_drop));
    Info("Imposed mass flow rate        =    " + to_string(out.mass_flow));
    Info("Imposed volumetric flow rate  =    " + to_string(out.volumetric_flow_rate));
  }
  Info("Inlet Velocity          =    " + to_string(out.inlet_velocity));
  Info("Outlet Velocity         =    " + to_string(out.outlet_velocity));
  Info("Inlet Pressure          =    " + to_string(cfg.atmospheric_pressure));
  Info("Outlet Pressure         =    " + to_string(cfg.atmospheric_pressure - out.pressure_drop));

  Info("==================== Stability and Accuracy Criteria =============");

  if (in.temporal_step == 1.0f && in.spatial_step == 1.0f && in.density == 1.0f)
  {
    Info("Initial macroscopic density, spatial step, and time step in the lattice framework are unity.");
  }
  else {
    Warning("Poor choice of spatial and temporal spacing, and density on the lattice. Should be non-dimensionalized to make Unity");
  }

  // Check stability and accuracy conditions

  if (in.knudsen_number > 0.1)
  {
    Warning("Knudsen numbers = " + to_string(in.knudsen_number) + " are not sufficiently small. Application of Lattice Boltzmann Equation is invalid. Redefine system parameters.");
  }
  else
  {
    Info("Knudsen numbers = " + to_string(in.knudsen_number) + " are sufficiently small to generate a valid hydrodynamic picture.");
  }

  in.reynolds_max = 4 * in.system_length / std::pow(in.speed_of_sound, 2);

  if (in.reynolds_number > in.reynolds_max)
  {
    Warning("Reynolds number = " + to_string(in.reynolds_number) + " is greater than maximum achievable Reynolds numbers. Reconsider simulation parameters");
  }
  else
  {
    Info("Reynolds number are less than maximum achievable Reynolds numbers (" + to_string(in.reynolds_max) + ") for this geometry.");
  };

  in.gamma = std::abs(std::log(cfg.temporal_step) / std::log(in.grid_spacing));
  if (in.gamma < 1)
  {
    Warning("Gamma     = " + to_string(in.gamma) + " < 1 : Compressibility errors are significant. Reconsider the temporal or spatial resolution.");
  }
  else
  {
    Info("Gamma     = " + to_string(in.gamma) + " > 1 : Compressibility errors are negligible");
  }

  if (in.mach_number > 0.3f)
  {
    Warning("Mach number = " + to_string(in.mach_number) + " is not small enough for weekly incompressible solution.");
  }
  else
  {
    Info("Mach number = " + to_string(in.mach_number) + " is small enough for weekly incompressible solution.");
  }

  in.essential_condition = 0.5f + 2.0f * (in.system_velocity) * in.reynolds_grid_number / std::pow(in.speed_of_sound, 2);
  if (in.relaxation_time < in.essential_condition)
  {
    Warning("Dimensionless relaxation parameter = " + to_string(in.relaxation_time) + " < stability threshold : increase viscosity or decrease grid spacing.");
  }
  else
  {
    Info("Dimensionless relaxation parameter = " + to_string(in.relaxation_time) + " is in the stable range.");
  }

  if (in.relaxation_time > in.essential_condition && in.relaxation_time < 2.0f && 2.0f * in.system_velocity <= 0.03)
  {
    Info("Relaxation time and maximum system velocity satisfy accuracy criterion");
  }
  else
  {
    Warning("Either relaxation parameter or the system velocity makes poor choice for an accurate solution.");
  }

  in.bc_accuracy = cfg.relaxation_time / cfg.temporal_step;
  Info("tau/Delta_t in physical units         = " + to_string(in.bc_accuracy));

  // Fill out the boundary conditions for simulation (should be non-dimensional)
  // Interior
  out.boundary_conditions[0] = cfg.atmospheric_pressure / conv_pressure;
  out.boundary_conditions[1] = cfg.initial_velocity / conv_velocity;
  out.boundary_conditions[2] = cfg.initial_velocity / conv_velocity;
  out.boundary_conditions[3] = cfg.initial_velocity / conv_velocity;
  // Inactive
  out.boundary_conditions[4] = cfg.atmospheric_pressure / conv_pressure;
  out.boundary_conditions[5] = 0;
  out.boundary_conditions[6] = 0;
  out.boundary_conditions[7] = 0;
  // Wall
  out.boundary_conditions[8] = cfg.atmospheric_pressure / conv_pressure;
  out.boundary_conditions[9] = 0;
  out.boundary_conditions[10] = 0;
  out.boundary_conditions[11] = 0;
  // Inflow
  out.boundary_conditions[12] = out.inlet_pressure / conv_pressure;
  out.boundary_conditions[13] = out.inlet_velocity / conv_velocity;
  out.boundary_conditions[14] = cfg.initial_velocity / conv_velocity;
  out.boundary_conditions[15] = cfg.initial_velocity / conv_velocity;
  // Outflow
  out.boundary_conditions[16] = out.outlet_pressure / conv_pressure;
  out.boundary_conditions[17] = out.outlet_velocity / conv_velocity;
  out.boundary_conditions[18] = cfg.initial_velocity / conv_velocity;
  out.boundary_conditions[19] = cfg.initial_velocity / conv_velocity;


  Info("Boundary Conditions in normalized units: (Pressure, X-Velocity, Y-Velocity, Z-Velocity)");
  Info("Interior:          " + to_string(out.boundary_conditions[0]) + "          " + to_string(out.boundary_conditions[1]) + "           " + to_string(out.boundary_conditions[2]) + "         " + to_string(out.boundary_conditions[3]));
  Info("Inactive:          " + to_string(out.boundary_conditions[4]) + "          " + to_string(out.boundary_conditions[5]) + "           " + to_string(out.boundary_conditions[6]) + "         " + to_string(out.boundary_conditions[7]));
  Info("Wall:              " + to_string(out.boundary_conditions[8]) + "          " + to_string(out.boundary_conditions[9]) + "           " + to_string(out.boundary_conditions[10]) + "         " + to_string(out.boundary_conditions[11]));
  Info("Inlet:             " + to_string(out.boundary_conditions[12]) + "          " + to_string(out.boundary_conditions[13]) + "           " + to_string(out.boundary_conditions[14]) + "         " + to_string(out.boundary_conditions[15]));
  Info("Outlet:            " + to_string(out.boundary_conditions[16]) + "          " + to_string(out.boundary_conditions[17]) + "           " + to_string(out.boundary_conditions[18]) + "         " + to_string(out.boundary_conditions[19]));


  // Setup the initial pressure and flow values for every cell
  int iCell = 0;
  for (size_t z = 0; z < in.dimensions[2]; z++)
  {
    for (size_t y = 0; y < in.dimensions[1]; y++)
    {
      for (size_t x = 0; x < in.dimensions[0]; x++)
      {
        int label = out.labels[iCell];
        if (label >= 0 && label < 99) // Wall and interior
        {
          out.temperature[iCell] = cfg.wall_temperature;                                          // Temperature
          out.pressure_flows_and_stress[iCell] = out.boundary_conditions[0];                      // Pressure
          out.pressure_flows_and_stress[out.num_cells + iCell] = out.boundary_conditions[1];      // X Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 2) + iCell] = out.boundary_conditions[2];// Y Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 3) + iCell] = out.boundary_conditions[3];// Z Flow Velocity
        }
        else if (label >= 101 && label < 200) // Inlet
        {
          out.temperature[iCell] = cfg.fluid_temperature;                                          // Temperature
          out.pressure_flows_and_stress[iCell] = out.boundary_conditions[12];                      // Pressure
          out.pressure_flows_and_stress[out.num_cells + iCell] = out.boundary_conditions[13];      // X Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 2) + iCell] = out.boundary_conditions[14];// Y Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 3) + iCell] = out.boundary_conditions[15];// Z Flow Velocity
        }
        else if (label >= 201 && label < 300) // Outlet
        {
          out.temperature[iCell] = cfg.wall_temperature;                                           // Temperature
          out.pressure_flows_and_stress[iCell] = out.boundary_conditions[16];                      // Pressure
          out.pressure_flows_and_stress[out.num_cells + iCell] = out.boundary_conditions[17];      // X Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 2) + iCell] = out.boundary_conditions[18];// Y Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 3) + iCell] = out.boundary_conditions[19];// Z Flow Velocity
        }
        else
        {
          out.temperature[iCell] = cfg.fluid_temperature;                                         // Temperature
          out.pressure_flows_and_stress[iCell] = out.boundary_conditions[4];                      // Pressure
          out.pressure_flows_and_stress[out.num_cells + iCell] = out.boundary_conditions[5];      // X Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 2) + iCell] = out.boundary_conditions[6];// Y Flow Velocity
          out.pressure_flows_and_stress[(out.num_cells * 3) + iCell] = out.boundary_conditions[7];// Z Flow Velocity
        }
        out.pressure_flows_and_stress[(out.num_cells * 4) + iCell] = 0.f;   //Sigma_xy
        out.pressure_flows_and_stress[(out.num_cells * 5) + iCell] = 0.f;   //Sigma_xz
        out.pressure_flows_and_stress[(out.num_cells * 6) + iCell] = 0.f;   //Sigma_yz
        out.pressure_flows_and_stress[(out.num_cells * 7) + iCell] = 0.f;   //Sigma_xx
        out.pressure_flows_and_stress[(out.num_cells * 8) + iCell] = 0.f;   //Sigma_yy
        out.pressure_flows_and_stress[(out.num_cells * 9) + iCell] = 0.f;   //Sigma_zz
        iCell++;
      }
    }
  }
  Info("LBM setup completed for "+cfg.name);
  // Only setup our structures and return
  if (cfg.setup_only)
    return true;

  in.pre_conditioning = false;

  cfg.print_stride = 100;    // Print error after these many iterations
  cfg.print_iteration = 0;   // Controls from which iteration to start printing

  if (in.pre_conditioning)
  {
    // Setup call to initiation microLBM call
    in.model_type = ModelType::Preconditioner;
    in.initialize_distribution_functions_from_cont_mech = true;
    in.return_distribution_functions = true; // Flag to return PDFs to host

    Info("Running initial microLBM as a pre-conditioner...");

    memset(out.previous_results, 0, out.num_cells * 3);// Reset our previous results

    for (size_t i = 0; i < out.num_cells * 3; i++)
    {
      out.previous_results[i] = out.pressure_flows_and_stress[num_cells + i];
    }

    status = microLBM(*this);
    if (status > 0)
    {
      Error("  Error occurred while running microLBM for pre-conditioner");
      return false;
    }
    else
    {
        for (int i = 0; i < out.nPDF; i++)
        {
            out.fPDF[i] *= 1/std::pow(lattice_speed,2);
        }
    }
  }

  // Advance time (run iterations of the LBM solver to converge to a solution)
  in.model_type = ModelType::Fluid;
  in.return_distribution_functions = true; // Flag to return PDFs to host
  cfg.print_iteration = 0;   // Controls from which iteration to start printing


  Info("Running microLBM on "+cfg.name+" for pressure and velocity...");

  memset(out.previous_results, 0, out.num_cells * 3);// Reset any previous velocity calculation
  for (size_t i = 0; i < out.num_cells * 3; i++)
  {
    out.previous_results[i] = out.pressure_flows_and_stress[num_cells + i];
  }

  timer.Start("LBM");
  // Set default value of PDF initialization flag based on whether the pre-conditioner was run or not
  in.initialize_distribution_functions_from_cont_mech = !in.pre_conditioning;
  status = microLBM(*this);

  if (status > 0)
  {
    Error("  Error occurred while running microLBM for pressure and velocity");
    return false;
  }
  out.compute_time_s = timer.GetElapsedTime_s("LBM");

  // Calculate temperature
  Info("Running microLBM for temperature...");
  for (size_t j = 0; j < out.num_cells; j++)
  {
    float p = out.pressure_flows_and_stress[j];
    out.temperature[out.num_cells + j] = out.pressure_flows_and_stress[out.num_cells + j] / p;        // X Flow Velocity / Pressure
    out.temperature[out.num_cells * 2 + j] = out.pressure_flows_and_stress[out.num_cells * 2 + j] / p;// Y Flow Velocity / Pressure
    out.temperature[out.num_cells * 3 + j] = out.pressure_flows_and_stress[out.num_cells * 3 + j] / p;// Z Flow Velocity / Pressure
  }

  iCell = 0;
  for (size_t z = 0; z < in.dimensions[2]; z++)
  {
    for (size_t y = 0; y < in.dimensions[1]; y++)
    {
      for (size_t x = 0; x < in.dimensions[0]; x++)
      {
        int label = out.labels[iCell];
        if (label >= 0 && label < 99) // Wall and interior
        {
          out.temperature[iCell] = cfg.wall_temperature;                                          // Temperature
        }
        else if (label >= 101 && label < 200) // Inlet
        {
          out.temperature[iCell] = cfg.fluid_temperature;                                          // Temperature
        }
        else if (label >= 201 && label < 300) // Outlet
        {
          out.temperature[iCell] = cfg.wall_temperature;                                           // Temperature
        }
        else
        {
          out.temperature[iCell] = cfg.fluid_temperature;                                         // Temperature
        }
        iCell++;
      }
    }
  }

  in.model_type = ModelType::Thermal;
  in.initialize_distribution_functions_from_cont_mech = true; // Flag whether to initialize PDFs from Q
  in.return_distribution_functions = true;        // Flag to return PDFs to host
  cfg.print_iteration = 0;   // Controls from which iteration to start printing


    memset(out.previous_results, 0, out.num_cells * 3);// Reset our previous results
  for (size_t i = 0; i < out.num_cells; i++)
  {
    out.previous_results[i] = out.temperature[i];
  }

  status = microLBM(*this);

  if (status > 0)
  {
    Error("  Error occurred while running microLBM for temperature");
    return false;
  }

  // Convert pressures, flows, stresses to the provided input units
  for (size_t i = 0; i < out.num_cells; i++)
  {
    float p = (out.pressure_flows_and_stress[i] * conv_pressure) - cfg.atmospheric_pressure;
    out.pressure_flows_and_stress[i] = p;
    out.pressure_flows_and_stress[out.num_cells + i] *= conv_velocity;
    out.pressure_flows_and_stress[out.num_cells * 2 + i] *= conv_velocity;
    out.pressure_flows_and_stress[out.num_cells * 3 + i] *= conv_velocity;
    out.pressure_flows_and_stress[out.num_cells * 4 + i] *= conv_pressure;
    out.pressure_flows_and_stress[out.num_cells * 5 + i] *= conv_pressure;
    out.pressure_flows_and_stress[out.num_cells * 6 + i] *= conv_pressure;
    out.pressure_flows_and_stress[out.num_cells * 7 + i] *= conv_pressure;
    out.pressure_flows_and_stress[out.num_cells * 8 + i] *= conv_pressure;
    out.pressure_flows_and_stress[out.num_cells * 9 + i] *= conv_pressure;
  }

  //  ===============================Compute Wall Shear Stress========================
  //  ================================================================================

  Info("================Computing Wall Shear Stress=============");
  int count = 0;  // count the wall sites where stress tensor is singular

  int   n_x, n_y, n_z;
  float n0_x, n0_y, n0_z;
  float n3_x, n3_y, n3_z;

  float Sxy, Sxz, Syz;
  float Sxx, Syy, Szz;
  float Sx, Sy, Sz;
  float Sx_new, Sy_new, Sz_new;

  float u_x, u_y, u_z;

  for (size_t i = 0; i < out.num_cells; ++i)
  {
    int label = out.labels[i];
    if (label > 0 && label <= 99) // wall boundary condition area
    {
      // -----------------------------------------------------------
      // Computing Shear Stress by Uing the Default Existing Normals
      //-------------------------------------------------------------
      Sxy = out.pressure_flows_and_stress[out.num_cells * 4 + i];
      Sxz = out.pressure_flows_and_stress[out.num_cells * 5 + i];
      Syz = out.pressure_flows_and_stress[out.num_cells * 6 + i];

      // Components of normals as obtained from the FindNormals function in LBM.cu
      n_x = VelocityDirections[label][2];
      n_y = VelocityDirections[label][1];
      n_z = VelocityDirections[label][0];

      //std::cout<<"label = "<<label<<" normal = "<<n_x<<","<<n_y<<","<<n_z<<"\n";

      Eigen::Vector3f Existing_Normal, Existing_Unit_Normal, Sigma_Default;
      Existing_Normal << float(n_x), float(n_y), float(n_z);
      Existing_Unit_Normal = Existing_Normal / Existing_Normal.norm();
      n0_x = Existing_Unit_Normal(0);
      n0_y = Existing_Unit_Normal(1);
      n0_z = Existing_Unit_Normal(2);

      float tmp = n0_x * n0_y * n0_z;

      Sx = Sxy * n0_y - Syz * tmp;
      Sy = Syz * n0_z - Sxz * tmp;
      Sz = Sxz * n0_x - Sxy * tmp;
      Sigma_Default << Sx, Sy, Sz;

      out.wall_shear_stress[i] = Sigma_Default.norm();  // This is the result using the existing normals

      // ---------------------------------------------------------------------
      // Computing Shear Stress by Using Wall Normals Algorithm by Stahl et. al.
      //-----------------------------------------------------------------------

      Sxx = out.pressure_flows_and_stress[out.num_cells * 7 + i];
      Syy = out.pressure_flows_and_stress[out.num_cells * 8 + i];
      Szz = out.pressure_flows_and_stress[out.num_cells * 9 + i];

      u_x = out.pressure_flows_and_stress[out.num_cells + i];
      u_y = out.pressure_flows_and_stress[out.num_cells * 2 + i];
      u_z = out.pressure_flows_and_stress[out.num_cells * 3 + i];

      Eigen::Matrix3f Stress_Tensor, An_Eigen_Base;
      Eigen::Vector3f Velocity, N1, N2, N3, New_Unit_Normal, An_Eigen_Set, Sigma_New;
      //Eigen::Vector3f norm_Eigen_Base;

      Stress_Tensor << Sxx, Sxy, Sxz,
        Sxy, Syy, Syz,
        Sxz, Syz, Szz;
      //std::cout<<Stress_Tensor<<"\n";
      // The if statement below checks if there are any wall sites for which analytical computation of the normal is possible.
      if (Stress_Tensor.determinant() >= 0.f && Stress_Tensor.determinant() <= 1E-20)
      {
        //std::cout<<"Sigma is singular\n";
        count++;
      }

      Velocity << u_x, u_y, u_z;
      N1 = Velocity / Velocity.norm();

      Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigensolver(Stress_Tensor);
      An_Eigen_Set = eigensolver.eigenvalues().cwiseAbs(); // Absolute of the Eigen values for Stress Tensor
      An_Eigen_Base = eigensolver.eigenvectors();          // Corresponding Unit Eigen Vectors of the Stress Tesnor (Bases for the Eigen Space)
     // norm_Eigen_Base <<  An_Eigen_Base.col(0).norm(),  An_Eigen_Base.col(1).norm(),  An_Eigen_Base.col(2).norm();  //Norm of Eigen vectors. Should be 1.

      // Finding the smallest Eigen value and its corresponding Eigen vector. This will be treated as N2.
      float Smallest_Eigen_Value = An_Eigen_Set.minCoeff();
      for (int j = 0; j < 3; j++)
      {
        if (An_Eigen_Set(j) == Smallest_Eigen_Value)
        {
          N2 = An_Eigen_Base.col(j);  // For being in the Bases of an Eigen Space, N2 is already a unit vector
        }
      }

      N3 = N1.cross(N2);
      New_Unit_Normal = N3 / N3.norm();

      n3_x = New_Unit_Normal(0);
      n3_y = New_Unit_Normal(1);
      n3_z = New_Unit_Normal(2);

      float tmp1 = n3_x * n3_y * n3_z;

      Sx_new = Sxy * n3_y - Syz * tmp1;
      Sy_new = Syz * n3_z - Sxz * tmp1;
      Sz_new = Sxz * n3_x - Sxy * tmp1;

      Sigma_New << Sx_new, Sy_new, Sz_new;

      // Final result using the newly computed normals
      out.wall_shear_stress[out.num_cells + i] = Sigma_New.norm();
    }
    else
    {
      out.wall_shear_stress[i] = 0.f;
      out.wall_shear_stress[out.num_cells + i] = 0.f;
    }
  }
  Info("Out of " + to_string(out.num_wall_cells) + " wall sites, viscous stress tensor is singular at =" + to_string(count) + " sites");
  Info("===============Finished Computing Wall Shear Stress=================");

  return true;
}

/* FindWalls: Given the array g with boundary codes:
     -1   = inactive node
      0   = active node
      >63 = flow condition node
   modify g such that any inactive node linked by a lattice direction to an active node becomes a wall
   node with the appropriate inward-pointing normal direction. Lattice direction numbering convention is
   encoded in least-significant 6 bits:
*/
bool LBM::Run::GetNormals()
{
  bool print_wall_sites = false;
  bool debug_wall_sites = false;

  int TRUE = 1, INTERIOR = 0, INACTIVE = -1; // FALSE=0
  int j[3], idx1[3], idx2[3], n, d;
  int NeighborInLattice1, NeighborInLattice2;
  size_t nloc1, loc;
  int nWallSites, nFlowSites, BCcategory;
  int FoundInteriorNeighbor, InactiveSite, BCFlowSiteNoNormal;
  int DefineWallSiteNormal, DefineFlowSiteNormal;                 // LBM model

  nWallSites = 0;
  nFlowSites = 0;
  if (debug_wall_sites)
  {
    Info("FindNormals:");
    Info("  L=" + to_string(in.dimensions));
  }
  for (j[0] = 0; j[0] < in.dimensions[0]; j[0]++)
  {
    for (j[1] = 0; j[1] < in.dimensions[1]; j[1]++)
    {
      for (j[2] = 0; j[2] < in.dimensions[2]; j[2]++)
      {
        loc = j[0] + in.dimensions[0] * (j[1] + in.dimensions[1] * j[2]);
        if (out.labels[loc] != INTERIOR)
        {
          // See if the inward-pointing normal direction must be defined at this site
          // Is this site inactive? (but perhaps a wall site for which a normal needs to be defined)
          InactiveSite = out.labels[loc] < 0;
          // Is this site on a flow condition boundary for which the normal has not yet been defined?
          BCFlowSiteNoNormal = out.labels[loc] % 100 == 0;
          BCcategory = out.labels[loc] / 100;
          if (debug_wall_sites)
          {
            Info("  Checking site at i=" + to_string(j[0]) + ", j=" + to_string(j[1]) + ", k=" + to_string(j[2]));
            Info("    bc=" + to_string(out.labels[loc]));
          }

          for (n = 1; n < nPDFmodels; n++)
          {
            // Find indices of neighboring lattice sites along each direction
            for (d = 0; d < 3; d++)
            {
              idx1[d] = j[d] + VelocityDirections[n][2 - d];
              idx2[d] = j[d] + 2 * VelocityDirections[n][2 - d];
            }
            NeighborInLattice1 = TRUE;
            NeighborInLattice2 = TRUE;
            for (d = 0; d < 3; d++)
            {
              NeighborInLattice1 = NeighborInLattice1 && (idx1[d] >= 0) && (idx1[d] < in.dimensions[d]);
              NeighborInLattice2 = NeighborInLattice2 && (idx2[d] >= 0) && (idx2[d] < in.dimensions[d]);
            }
            if (debug_wall_sites)
            {
              Info("  Along direction nr " + to_string(n) + " idx= " + to_string(idx1[0]) + to_string(idx1[1]) + to_string(idx1[2]));
              Info("    NeighborInLattice1=" + to_string(NeighborInLattice1));
            }
            if (NeighborInLattice1)
            {
              // idx1[:] defines a valid interior site
              nloc1 = idx1[0] + in.dimensions[0] * (idx1[1] + in.dimensions[1] * idx1[2]);
              // Is this neighbor lattice site in the interior?
              FoundInteriorNeighbor = out.labels[nloc1] == INTERIOR;
              // Is this a wall site, i.e., an inactive site with an interior neighbor?
              // If so, an interior pointing normal direction must be defined
              DefineWallSiteNormal = InactiveSite && FoundInteriorNeighbor;
              // Does this flow BC site have an interior neighbor?
              // If so, an interior pointing normal direction must be defined
              DefineFlowSiteNormal = BCFlowSiteNoNormal && FoundInteriorNeighbor;
              if (debug_wall_sites)
              {
                Info("  FoundInteriorNeighbor=" + to_string(FoundInteriorNeighbor));
                Info("  DefineWallSiteNormal=" + to_string(DefineWallSiteNormal));
                Info("  DefineFlowSiteNormal=" + to_string(DefineFlowSiteNormal));
              }
              if (DefineWallSiteNormal)
              {
                out.labels[loc] = n;
                nWallSites++;
                if (print_wall_sites)
                {
                  Info("  Found wall node at i=" + to_string(j[0]) + ", j=" + to_string(j[1]) + ", k=" + to_string(j[2]));
                  Info("  n=" + to_string(n) + " normal=" + to_string(VelocityDirections[n]) + " G=" + to_string(out.labels[loc]));
                  Info("  Neighbor at i=" + to_string(idx1[0]) + ", j=" + to_string(idx1[1]) + ", k=" + to_string(idx1[2]));
                  Info("    G=" + to_string(out.labels[nloc1]));
                }
                break;
              }
              if (DefineFlowSiteNormal)
              {
                // Store lattice direction corresponding to interior normal
                out.labels[loc] = BCcategory * 100 + n;
                nFlowSites++;
                if (!NeighborInLattice2)
                {
                  Error("  Grid does not contain two nodes interior to flow boundary condition");
                  exit(-2);
                }
                if (print_wall_sites)
                {
                  if (BCcategory == 1)
                  {
                    Info("  Have defined inflow normal at i=" + to_string(j[0]) + ", j=" + to_string(j[1]) + ", k=" + to_string(j[2]));
                    Info("    n=" + to_string(n) + " normal=" + to_string(VelocityDirections[n]) + " G=" + to_string(out.labels[loc]));
                  }
                  else
                  {
                    Info("  Have defined outflow normal at i=" + to_string(j[0]) + ", j=" + to_string(j[1]) + ", k=" + to_string(j[2]));
                    Info("    n=" + to_string(n) + " normal=" + to_string(VelocityDirections[n]) + " G=" + to_string(out.labels[loc]));
                  }
                }
                break;
              }
            }
          }
          // Flow BC sites not identified as connected to an interior site are set inactive
          if ((out.labels[loc] >= 100) && (out.labels[loc] % 100 == 0))
            out.labels[loc] = INACTIVE;
        }
      }
    }
  }
  Info("FindNormals has defined normals at " + to_string(nWallSites) + " wall sites, " + to_string(nFlowSites) + " flow sites");

  return true;
}

void LBM::Configuration::ToString(std::ostream& str) const
{
  str << "LBM Configuration : \n";
  str << "  atmospheric_pressure = " << atmospheric_pressure << "\n";
  str << "  density = " << density << "\n";
  str << "  fluid_temperature = " << fluid_temperature << "\n";
  str << "  initial_velocity = " << initial_velocity << "\n";
  str << "  speed_of_sound = " << speed_of_sound << "\n";
  str << "  viscosity = " << viscosity << "\n";
  str << "  dimensionless viscosity = " << unitless_viscosity << "\n";
  str << "  initial flow rate = " << target_flow_rate << "\n";
  str << "  thermal_coefficient = " << thermal_coefficient << "\n";
  str << "  wall_temperature = " << wall_temperature << "\n";
  str << "  inlet_boundary_condition = " << (inlet_boundary_condition == BoundaryCondition::Pressure ? "Pressure" : "Volumetric Flow Rate") << "\n";
  str << "  inlet_pressure = " << inlet_pressure;
  str << "  inlet_volumetric_flow_rate = " << inlet_volumetric_flow_rate;
  str << "  outlet_boundary_condition = " << (outlet_boundary_condition == BoundaryCondition::Pressure ? "Pressure" : "Volumetric Flow Rate") << "\n";
  str << "  outlet_pressure = " << outlet_pressure;
  str << "  outlet_volumetric_flow_rate = " << outlet_volumetric_flow_rate;
  str << "  max_iterations = " << max_iterations;
  str << "  tolerance = " << tolerance;
}

void LBM::FlowMetrics::ToString(std::ostream& str) const
{
  str << "Flow Metrics : \n";
  str << "  Flow Rate = " << flow_rate << "\n";
  str << "  Resistance = " << resistance << "\n";
  str << "  Heat Flux = " << heat_flux << "\n";
  str << "  wallShear = " << wall_shear;
}
