From e1ebbce812550a0711f92c351f52930ef17ba2d6 Mon Sep 17 00:00:00 2001
From: Zhe Rex Li <zhe.rex.li@gmail.com>
Date: Tue, 12 Apr 2022 11:40:39 +1000
Subject: [PATCH] initialize nernst-planck simulator; to be built subject to
 debugging

---
 models/MultiPhysController.cpp         |  46 +++++++++
 models/MultiPhysController.h           |   2 +-
 tests/CMakeLists.txt                   |   1 +
 tests/lbpm_nernst_planck_simulator.cpp | 135 +++++++++++++++++++++++++
 4 files changed, 183 insertions(+), 1 deletion(-)
 create mode 100644 tests/lbpm_nernst_planck_simulator.cpp

diff --git a/models/MultiPhysController.cpp b/models/MultiPhysController.cpp
index 8935bca7..9120471b 100644
--- a/models/MultiPhysController.cpp
+++ b/models/MultiPhysController.cpp
@@ -151,6 +151,52 @@ vector<int> ScaLBL_Multiphys_Controller::getIonNumIter_PNP_coupling(
     return num_iter_ion;
 }
 
+vector<int> ScaLBL_Multiphys_Controller::getIonNumIter_NernstPlanck_coupling(
+    const vector<double> &IonTimeConv) {
+    //Return number of internal iterations for the Ion transport solver
+    vector<int> num_iter_ion;
+    vector<double>::iterator it_max =
+        max_element(IonTimeConv.begin(), IonTimeConv.end());
+    unsigned int idx_max = distance(IonTimeConv.begin(), it_max);
+    if (idx_max == 0) {
+        num_iter_ion.push_back(2);
+        for (unsigned int idx = 1; idx < IonTimeConv.size(); idx++) {
+            double temp =
+                2 * TimeConv[idx_max] /
+                TimeConv
+                    [idx]; //the factor 2 is the number of iterations for the element has max time_conv
+            num_iter_ion.push_back(int(round(temp / 2) * 2));
+        }
+    } else if (idx_max == IonTimeConv.size() - 1) {
+        for (unsigned int idx = 0; idx < TimeConv.size() - 1; idx++) {
+            double temp =
+                2 * TimeConv[idx_max] /
+                TimeConv
+                    [idx]; //the factor 2 is the number of iterations for the element has max time_conv
+            num_iter_ion.push_back(int(round(temp / 2) * 2));
+        }
+        num_iter_ion.push_back(2);
+    } else {
+        for (unsigned int idx = 0; idx < idx_max; idx++) {
+            double temp =
+                2 * TimeConv[idx_max] /
+                TimeConv
+                    [idx]; //the factor 2 is the number of iterations for the element has max time_conv
+            num_iter_ion.push_back(int(round(temp / 2) * 2));
+        }
+        num_iter_ion.push_back(2);
+        for (unsigned int idx = idx_max + 1; idx < IonTimeConv.size(); idx++) {
+            double temp =
+                2 * TimeConv[idx_max] /
+                TimeConv
+                    [idx]; //the factor 2 is the number of iterations for the element has max time_conv
+            num_iter_ion.push_back(int(round(temp / 2) * 2));
+        }
+    }
+    return num_iter_ion;
+}
+
+
 void ScaLBL_Multiphys_Controller::getTimeConvMax_PNP_coupling(
     double StokesTimeConv, const vector<double> &IonTimeConv) {
     //Return maximum of the time converting factor from Stokes and ion solvers
diff --git a/models/MultiPhysController.h b/models/MultiPhysController.h
index 13aa05a2..d41a6db8 100644
--- a/models/MultiPhysController.h
+++ b/models/MultiPhysController.h
@@ -28,7 +28,7 @@ public:
                                       const vector<double> &IonTimeConv);
     vector<int> getIonNumIter_PNP_coupling(double StokesTimeConv,
                                            const vector<double> &IonTimeConv);
-    //void getIonNumIter_PNP_coupling(double StokesTimeConv,vector<double> &IonTimeConv,vector<int> &IonTimeMax);
+    vector<int> getIonNumIter_NernstPlanck_coupling(const vector<double> &IonTimeConv);
     void getTimeConvMax_PNP_coupling(double StokesTimeConv,
                                      const vector<double> &IonTimeConv);
 
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 5e4441e3..2ff1aa88 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -6,6 +6,7 @@ ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_electrokinetic_SingleFluid_simulator )
+ADD_LBPM_EXECUTABLE( lbpm_nernst_planck_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_freelee_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_freelee_SingleFluidBGK_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
diff --git a/tests/lbpm_nernst_planck_simulator.cpp b/tests/lbpm_nernst_planck_simulator.cpp
new file mode 100644
index 00000000..bc697bd9
--- /dev/null
+++ b/tests/lbpm_nernst_planck_simulator.cpp
@@ -0,0 +1,135 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/stat.h>
+#include <iostream>
+#include <exception>
+#include <stdexcept>
+#include <fstream>
+#include <math.h>
+
+#include "models/IonModel.h"
+#include "models/PoissonSolver.h"
+#include "models/MultiPhysController.h"
+#include "common/Utilities.h"
+#include "analysis/ElectroChemistry.h"
+
+using namespace std;
+
+int main(int argc, char **argv)
+{
+    // Initialize MPI and error handlers
+    	Utilities::startup( argc, argv );
+    	Utilities::MPI comm( MPI_COMM_WORLD );
+    	int rank = comm.getRank();
+    	int nprocs = comm.getSize();
+
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+        if (rank == 0){
+            printf("*************************************************\n");
+            printf("Running LBPM Nernst-Planck solver \n");
+            printf("*************************************************\n");
+        }
+    	// Initialize compute device
+    	int device=ScaLBL_SetDevice(rank);
+        NULL_USE( device );
+    	ScaLBL_DeviceBarrier();
+    	comm.barrier();
+
+        PROFILE_ENABLE(1);
+        //PROFILE_ENABLE_TRACE();
+        //PROFILE_ENABLE_MEMORY();
+        PROFILE_SYNCHRONIZE();
+        PROFILE_START("Main");
+        Utilities::setErrorHandlers();
+
+        auto filename = argv[1];
+        ScaLBL_IonModel IonModel(rank,nprocs,comm);
+        ScaLBL_Poisson PoissonSolver(rank,nprocs,comm); 
+        ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling
+        
+        // Load controller information
+        Study.ReadParams(filename);
+
+        // Load user input database files for Ion solvers
+        IonModel.ReadParams(filename);
+        IonModel.SetDomain();    
+        IonModel.ReadInput();    
+        IonModel.Create();      
+        
+        // Create analysis object
+        //ElectroChemistryAnalyzer Analysis(IonModel.Dm);
+
+
+        //IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
+        IonModel.timestepMax = Study.getIonNumIter_NernstPlanck_coupling(IonModel.time_conv);
+        IonModel.Initialize();   
+        
+        // Get maximal time converting factor based on Sotkes and Ion solvers
+        //Study.getTimeConvMax_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
+        // Get time conversion factor for the main iteration loop in electrokinetic single fluid simulator
+        Study.time_conv_MainLoop = IonModel.timestepMax[0]*IonModel.time_conv[0];
+
+        //----------------------------------- print out for debugging ------------------------------------------//
+        if (rank==0){
+            for (size_t i=0;i<IonModel.timestepMax.size(),i++){
+                printf("Main loop time_conv computed from %ith ion: %.5g[s/lt]"%(IonModel.timestepMax[i]*IonModel.time_conv[i]));
+            }
+        }
+        //----------------------------------- print out for debugging ------------------------------------------//
+
+        // Initialize LB-Poisson model
+        PoissonSolver.ReadParams(filename);
+        PoissonSolver.SetDomain();    
+        PoissonSolver.ReadInput();    
+        PoissonSolver.Create();       
+        PoissonSolver.Initialize(Study.time_conv_MainLoop);   
+
+
+        if (rank == 0){
+            printf("********************************************************\n");
+            printf("Key Summary of LBPM electrokinetic single-fluid solver \n");
+            printf("   1. Max LB Timestep: %i [lt]\n", Study.timestepMax);
+            printf("   2. Time conversion factor per LB Timestep: %.6g [sec/lt]\n",Study.time_conv_MainLoop);
+            printf("   3. Max Physical Time: %.6g [sec]\n",Study.timestepMax*Study.time_conv_MainLoop);
+            printf("********************************************************\n");
+        }
+
+        int timestep=0;
+        while (timestep < Study.timestepMax){
+            
+            timestep++;
+            PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
+            //StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
+            IonModel.Run(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField); //solve for ion transport and electric potential
+            
+
+            //if (timestep%Study.analysis_interval==0){
+            //	Analysis.Basic(IonModel,PoissonSolver,StokesModel,timestep);
+            //}
+            if (timestep%Study.visualization_interval==0){
+                //Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep);
+            	PoissonSolver.getElectricPotential_debug(timestep);
+                PoissonSolver.getElectricField_debug(timestep);
+                IonModel.getIonConcentration_debug(timestep);
+                //StokesModel.getVelocity(timestep);
+            	
+            }
+        }
+
+        //if (rank==0) printf("Save simulation raw data at maximum timestep\n");
+    	//Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep);
+
+        if (rank==0) printf("Maximum LB timestep = %i is reached and the simulation is completed\n",Study.timestepMax);
+        if (rank==0) printf("*************************************************************\n");
+
+        PROFILE_STOP("Main");
+        PROFILE_SAVE("lbpm_nernst_planck_simulator",1);
+        // ****************************************************
+
+    } // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+  Utilities::shutdown();
+}
+
+