From 613a43c309a288c02a9b5ffdff3c9f7f434df616 Mon Sep 17 00:00:00 2001
From: James McClure <jemcclur@gmail.com>
Date: Wed, 26 Oct 2022 16:17:23 -0400
Subject: [PATCH] vis capabilities for poisson, default d3q19

---
 models/IonModel.cpp      |  1 -
 models/PoissonSolver.cpp | 50 +++++++++++++++++++++++++++++++++++++++-
 models/PoissonSolver.h   |  2 ++
 3 files changed, 51 insertions(+), 2 deletions(-)

diff --git a/models/IonModel.cpp b/models/IonModel.cpp
index 21f55225..4e389de3 100644
--- a/models/IonModel.cpp
+++ b/models/IonModel.cpp
@@ -1549,7 +1549,6 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
 
 			IonMembrane->IonTransport(&fq[ic * Np * 7],&Ci[ic * Np]);
 
-
 			/*           if (BoundaryConditionSolid == 1) {
                 //TODO IonSolid may also be species-dependent
                 ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic * Np * 7], IonSolid);
diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp
index 4f611204..89bfb985 100644
--- a/models/PoissonSolver.cpp
+++ b/models/PoissonSolver.cpp
@@ -86,7 +86,7 @@ void ScaLBL_Poisson::ReadParams(string filename){
 	}
     //'tolerance_method' can be {"MSE","MSE_max"}
 	tolerance_method = electric_db->getWithDefault<std::string>( "tolerance_method", "MSE" );
-	lattice_scheme = electric_db->getWithDefault<std::string>( "lattice_scheme", "D3Q7" );
+	lattice_scheme = electric_db->getWithDefault<std::string>( "lattice_scheme", "D3Q19" );
 	if (electric_db->keyExists( "epsilonR" )){
 		epsilonR = electric_db->getScalar<double>( "epsilonR" );
 	}
@@ -1320,6 +1320,54 @@ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
     }
 }
 
+void ScaLBL_Poisson::WriteVis( int timestep) {
+
+    auto vis_db = db->getDatabase("Visualization");
+    char VisName[40];
+    auto format = vis_db->getWithDefault<string>( "format", "hdf5" );
+
+    DoubleArray ElectricalPotential(Nx, Ny, Nz);
+    std::vector<IO::MeshDataStruct> visData;
+    fillHalo<double> fillData(Dm->Comm, Dm->rank_info,
+                              {Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2}, {1, 1, 1},
+                              0, 1);
+
+    IO::initialize("",format,"false");
+    // Create the MeshDataStruct    
+    visData.resize(1);
+
+    visData[0].meshName = "domain";
+    visData[0].mesh =
+        std::make_shared<IO::DomainMesh>(Dm->rank_info, Dm->Nx - 2, Dm->Ny - 2,
+                                         Dm->Nz - 2, Dm->Lx, Dm->Ly, Dm->Lz);
+    //electric potential
+    auto ElectricPotentialVar = std::make_shared<IO::Variable>();
+
+    //--------------------------------------------------------------------------------------------------------------------
+
+    //-------------------------------------Create Names for Variables------------------------------------------------------
+    if (vis_db->getWithDefault<bool>("save_electric_potential", true)) {
+        ElectricPotentialVar->name = "ElectricPotential";
+        ElectricPotentialVar->type = IO::VariableType::VolumeVariable;
+        ElectricPotentialVar->dim = 1;
+        ElectricPotentialVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2);
+        visData[0].vars.push_back(ElectricPotentialVar);
+    }
+    //--------------------------------------------------------------------------------------------------------------------
+
+    //------------------------------------Save All Variables--------------------------------------------------------------
+    if (vis_db->getWithDefault<bool>("save_electric_potential", true)) {
+        ASSERT(visData[0].vars[0]->name == "ElectricPotential");
+        getElectricPotential(ElectricalPotential);
+        Array<double> &ElectricPotentialData = visData[0].vars[0]->data;
+        fillData.copy(ElectricalPotential, ElectricPotentialData);
+    }
+
+    if (vis_db->getWithDefault<bool>("write_silo", true))
+        IO::writeData(timestep, visData, Dm->Comm);
+    //--------------------------------------------------------------------------------------------------------------------
+}
+
 //void ScaLBL_Poisson::SolveElectricField(){
 //    ScaLBL_Comm_Regular->SendHalo(Psi);
 //    ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid,
diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h
index faee8c01..c6672cec 100644
--- a/models/PoissonSolver.h
+++ b/models/PoissonSolver.h
@@ -22,6 +22,7 @@
 #ifndef ScaLBL_POISSON_INC
 #define ScaLBL_POISSON_INC
 
+
 class ScaLBL_Poisson {
 public:
     ScaLBL_Poisson(int RANK, int NP, const Utilities::MPI &COMM);
@@ -43,6 +44,7 @@ public:
                           DoubleArray &Values_z);
     void getElectricField_debug(int timestep);
     void Checkpoint();
+    void WriteVis( int timestep);
 
     void DummyChargeDensity(); //for debugging