OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

CPU only;finish preliminary work on testing Poisson and Ion models and their coupling

Browse Source
This commit is contained in:
Rex Zhe Li
2020-09-11 22:56:00 -04:00
parent 20c8cc9c3b
commit 11e9af54b6
18 changed files with 2453 additions and 226 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
common/ScaLBL.cpp
View File

@@ -2086,3 +2086,25 @@ void ScaLBL_Communicator::Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout){
ScaLBL_Poisson_D3Q7_BC_Z(dvcSendList_Z, Map, Psi, Vout, sendCount_Z);
}
}
void ScaLBL_Communicator::D3Q7_Ion_Concentration_BC_z(int *neighborList, double *fq, double Cin, int time){
if (kproc == 0) {
if (time%2==0){
ScaLBL_D3Q7_AAeven_Ion_Concentration_BC_z(dvcSendList_z, fq, Cin, sendCount_z, N);
}
else{
ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_z(neighborList, dvcSendList_z, fq, Cin, sendCount_z, N);
}
}
}
void ScaLBL_Communicator::D3Q7_Ion_Concentration_BC_Z(int *neighborList, double *fq, double Cout, int time){
if (kproc == nprocz-1){
if (time%2==0){
ScaLBL_D3Q7_AAeven_Ion_Concentration_BC_Z(dvcSendList_Z, fq, Cout, sendCount_Z, N);
}
else{
ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_Z(neighborList, dvcSendList_Z, fq, Cout, sendCount_Z, N);
}
}
}

16
common/ScaLBL.h
View File

@@ -79,17 +79,15 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, double zi, double rlx, double Vt, int start, int finish, int Np);
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, double zi, double rlx, double Vt, int start, int finish, int Np);
double Di, int zi, double rlx, double Vt, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);
extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np);
extern "C" void ScaLBL_IonConcentration_Phys(double *Den, double h, int ion_component, int start, int finish, int Np);
// LBM Poisson solver
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
@@ -206,6 +204,14 @@ extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi, doubl
extern "C" void ScaLBL_Poisson_D3Q7_BC_Z(int *list, int *Map, double *Psi, double Vout, int count);
extern "C" void ScaLBL_D3Q7_AAeven_Ion_Concentration_BC_z(int *list, double *dist, double Cin, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAeven_Ion_Concentration_BC_Z(int *list, double *dist, double Cout, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_z(int *d_neighborList, int *list, double *dist, double Cin, int count, int Np);
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_Z(int *d_neighborList, int *list, double *dist, double Cout, int count, int Np);
class ScaLBL_Communicator{
public:
//......................................................................................
@@ -269,6 +275,8 @@ public:
void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time);
void Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin);
void Poisson_D3Q7_BC_Z(int *Map, double *Psi, double Vout);
void D3Q7_Ion_Concentration_BC_z(int *neighborList, double *fq, double Cin, int time);
void D3Q7_Ion_Concentration_BC_Z(int *neighborList, double *fq, double Cout, int time);
// Debugging and unit testing functions
void PrintD3Q19();

86
cpu/D3Q7BC.cpp
View File

@@ -137,3 +137,89 @@ extern "C" void ScaLBL_Poisson_D3Q7_BC_Z(int *list, int *Map, double *Psi, doubl
Psi[nm] = Vout;
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion_Concentration_BC_z(int *list, double *dist, double Cin, int count, int Np){
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
double f1 = dist[2*Np+n];
double f2 = dist[1*Np+n];
double f3 = dist[4*Np+n];
double f4 = dist[3*Np+n];
double f6 = dist[5*Np+n];
//...................................................
double f5 = Cin - (f0+f1+f2+f3+f4+f6);
dist[6*Np+n] = f5;
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion_Concentration_BC_Z(int *list, double *dist, double Cout, int count, int Np){
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
double f1 = dist[2*Np+n];
double f2 = dist[1*Np+n];
double f3 = dist[4*Np+n];
double f4 = dist[3*Np+n];
double f5 = dist[6*Np+n];
//...................................................
double f6 = Cout - (f0+f1+f2+f3+f4+f5);
dist[5*Np+n] = f6;
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_z(int *d_neighborList, int *list, double *dist, double Cin, int count, int Np){
int nread,nr5;
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
nread = d_neighborList[n];
double f1 = dist[nread];
nread = d_neighborList[n+2*Np];
double f3 = dist[nread];
nread = d_neighborList[n+Np];
double f2 = dist[nread];
nread = d_neighborList[n+3*Np];
double f4 = dist[nread];
nread = d_neighborList[n+5*Np];
double f6 = dist[nread];
// Unknown distributions
nr5 = d_neighborList[n+4*Np];
double f5 = Cin - (f0+f1+f2+f3+f4+f6);
dist[nr5] = f5;
}
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion_Concentration_BC_Z(int *d_neighborList, int *list, double *dist, double Cout, int count, int Np){
int nread,nr6;
for (int idx=0; idx<count; idx++){
int n = list[idx];
double f0 = dist[n];
nread = d_neighborList[n];
double f1 = dist[nread];
nread = d_neighborList[n+2*Np];
double f3 = dist[nread];
nread = d_neighborList[n+4*Np];
double f5 = dist[nread];
nread = d_neighborList[n+Np];
double f2 = dist[nread];
nread = d_neighborList[n+3*Np];
double f4 = dist[nread];
// unknown distributions
nr6 = d_neighborList[n+5*Np];
double f6 = Cout - (f0+f1+f2+f3+f4+f5);
dist[nr6] = f6;
}
}

56
cpu/Ion.cpp
View File

@@ -81,7 +81,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_IonConcentration(double *dist, double *Den, i
}
extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, double zi, double rlx, double Vt, int start, int finish, int Np){
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
@@ -126,31 +126,31 @@ extern "C" void ScaLBL_D3Q7_AAodd_Ion(int *neighborList, double *dist, double *D
f6 = dist[nr6];
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*Ci;
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
// q = 1
dist[nr2] = f1*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(ux+uEPx));
dist[nr2] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx));
// q=2
dist[nr1] = f2*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(ux+uEPx));
dist[nr1] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx));
// q = 3
dist[nr4] = f3*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uy+uEPy));
dist[nr4] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy));
// q = 4
dist[nr3] = f4*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uy+uEPy));
dist[nr3] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy));
// q = 5
dist[nr6] = f5*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uz+uEPz));
dist[nr6] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz));
// q = 6
dist[nr5] = f6*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uz+uEPz));
dist[nr5] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz));
}
}
extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Velocity, double *ElectricField,
double Di, double zi, double rlx, double Vt, int start, int finish, int Np){
double Di, int zi, double rlx, double Vt, int start, int finish, int Np){
int n;
double Ci;
double ux,uy,uz;
@@ -181,25 +181,25 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *Veloci
f6 = dist[5*Np+n];
// q=0
dist[n] = f0*(1.0-rlx)+rlx*0.3333333333333333*Ci;
dist[n] = f0*(1.0-rlx)+rlx*0.25*Ci;
// q = 1
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(ux+uEPx));
dist[1*Np+n] = f1*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(ux+uEPx));
// q=2
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(ux+uEPx));
dist[2*Np+n] = f2*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(ux+uEPx));
// q = 3
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uy+uEPy));
dist[3*Np+n] = f3*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uy+uEPy));
// q = 4
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uy+uEPy));
dist[4*Np+n] = f4*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uy+uEPy));
// q = 5
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0+4.5*(uz+uEPz));
dist[5*Np+n] = f5*(1.0-rlx) + rlx*0.125*Ci*(1.0+4.0*(uz+uEPz));
// q = 6
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.1111111111111111*Ci*(1.0-4.5*(uz+uEPz));
dist[6*Np+n] = f6*(1.0-rlx) + rlx*0.125*Ci*(1.0-4.0*(uz+uEPz));
}
@@ -209,13 +209,13 @@ extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit,
{
int n;
for (n=0; n<Np; n++){
dist[0*Np+n] = 0.3333333333333333*DenInit;
dist[1*Np+n] = 0.1111111111111111*DenInit;
dist[2*Np+n] = 0.1111111111111111*DenInit;
dist[3*Np+n] = 0.1111111111111111*DenInit;
dist[4*Np+n] = 0.1111111111111111*DenInit;
dist[5*Np+n] = 0.1111111111111111*DenInit;
dist[6*Np+n] = 0.1111111111111111*DenInit;
dist[0*Np+n] = 0.25*DenInit;
dist[1*Np+n] = 0.125*DenInit;
dist[2*Np+n] = 0.125*DenInit;
dist[3*Np+n] = 0.125*DenInit;
dist[4*Np+n] = 0.125*DenInit;
dist[5*Np+n] = 0.125*DenInit;
dist[6*Np+n] = 0.125*DenInit;
Den[n] = DenInit;
}
}
@@ -236,13 +236,3 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
}
}
extern "C" void ScaLBL_IonConcentration_Phys(double *Den, double h, int ion_component, int start, int finish, int Np){
//h: resolution [um/lu]
int n;
double Ci;
for (n=start; n<finish; n++){
Ci = Den[n+ion_component*Np];
Den[n+ion_component*Np] = Ci/(h*h*h*1.0e-18);
}
}

1203
gpu/D3Q7BC.cu Normal file
View File

File diff suppressed because it is too large Load Diff

557
models/IonModel.cpp
View File

@@ -1,6 +1,7 @@
/*
* Dilute Ion Transport LBM Model
*/
#include <algorithm>
#include "models/IonModel.h"
#include "analysis/distance.h"
#include "common/ReadMicroCT.h"
@@ -8,6 +9,7 @@
ScaLBL_IonModel::ScaLBL_IonModel(int RANK, int NP, MPI_Comm COMM):
rank(RANK),nprocs(NP),timestep(0),timestepMax(0),time_conv(0),kb(0),electron_charge(0),T(0),Vt(0),k2_inv(0),h(0),
tolerance(0),number_ion_species(0),Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),
fluidVelx_dummy(0),fluidVely_dummy(0),fluidVelz_dummy(0),
BoundaryCondition(0),BoundaryConditionSolid(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
@@ -16,19 +18,13 @@ ScaLBL_IonModel::~ScaLBL_IonModel(){
}
void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stokes,double time_conv_Stokes){
void ScaLBL_IonModel::ReadParams(string filename,vector<int> &num_iter){
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
ion_db = db->getDatabase( "Ions" );
//------ Load number of iteration from multiphysics controller ------//
timestepMax = num_iter;
//compute time conversion factor for ion model
time_conv = num_iter_Stokes*time_conv_Stokes/num_iter;
//-------------------------------------------------------------------//
// Universal constant
kb = 1.38e-23;//Boltzmann constant;unit [J/K]
electron_charge = 1.6e-19;//electron charge;unit [C]
@@ -36,26 +32,30 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
//---------------------- Default model parameters --------------------------//
T = 300.0;//temperature; unit [K]
Vt = kb*T/electron_charge;//thermal voltage; unit [Vy]
k2_inv = 4.5;//speed of sound for D3Q7 lattice
k2_inv = 4.0;//speed of sound for D3Q7 lattice
h = 1.0;//resolution; unit: um/lu
tolerance = 1.0e-8;
number_ion_species = 1;
tau.push_back(1.0);
IonDiffusivity.push_back(1.0e-9);//user-input diffusivity has physical unit [m^2/sec]
IonValence.push_back(1);//algebraic valence charge
IonConcentration.push_back(1.0e-3);//user-input ion concentration has physical unit [mol/m^3]
//deltaT.push_back(1.0);
//tau.push_back(0.5+k2_inv*deltaT[0]*IonDiffusivity[0]);
tau.push_back(0.5+k2_inv*time_conv/(h*1.0e-6)/(h*1.0e-6)*IonDiffusivity[0]);
Cin.push_back(1.0e-3);//user-input inlet boundary ion concentration;unit [mol/m^3]
Cout.push_back(1.0e-3);//user-input outlet boundary ion concentration;unit [mol/m^3]
//tau.push_back(0.5+k2_inv*time_conv/(h*1.0e-6)/(h*1.0e-6)*IonDiffusivity[0]);
time_conv.push_back((tau[0]-0.5)/k2_inv*(h*h*1.0e-12)/IonDiffusivity[0]);
fluidVelx_dummy = 0.0;//for debugging, unit [m/sec]
fluidVely_dummy = 0.0;//for debugging, unit [m/sec]
fluidVelz_dummy = 0.0;//for debugging, unit [m/sec]
Ex_dummy = 0.0;//for debugging, unit [V/m]
Ey_dummy = 0.0;//for debugging, unit [V/m]
Ez_dummy = 0.0;//for debugging, unit [V/m]
//--------------------------------------------------------------------------//
// Read domain parameters
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
}
BoundaryCondition = 0;
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
// LB-Ion Model parameters
//if (ion_db->keyExists( "timestepMax" )){
@@ -69,29 +69,63 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
//re-calculate thermal voltage
Vt = kb*T/electron_charge;//thermal voltage; unit [Vy]
}
if (ion_db->keyExists( "FluidVelDummy" )){
fluidVelx_dummy = ion_db->getVector<double>( "FluidVelDummy" )[0];
fluidVely_dummy = ion_db->getVector<double>( "FluidVelDummy" )[1];
fluidVelz_dummy = ion_db->getVector<double>( "FluidVelDummy" )[2];
}
if (ion_db->keyExists( "ElectricFieldDummy" )){
Ex_dummy = ion_db->getVector<double>( "ElectricFieldDummy" )[0];
Ey_dummy = ion_db->getVector<double>( "ElectricFieldDummy" )[1];
Ez_dummy = ion_db->getVector<double>( "ElectricFieldDummy" )[2];
}
if (ion_db->keyExists( "number_ion_species" )){
number_ion_species = ion_db->getScalar<int>( "number_ion_species" );
}
//------ Load number of iteration from multiphysics controller ------//
if (num_iter.size()!=number_ion_species){
ERROR("Error: number_ion_species and num_iter_Ion_List (from Multiphysics) must be of the same length! \n");
}
else{
timestepMax.assign(num_iter.begin(),num_iter.end());
}
//-------------------------------------------------------------------//
if (ion_db->keyExists("tauList")){
tau.clear();
tau = ion_db->getVector<double>( "tauList" );
vector<double>Di = ion_db->getVector<double>( "IonDiffusivityList" );//temp storing ion diffusivity in physical unit
if (tau.size()!=number_ion_species || Di.size()!=number_ion_species){
ERROR("Error: number_ion_species, tauList and IonDiffusivityList must be of the same length! \n");
}
else{
time_conv.clear();
for (int i=0; i<tau.size();i++){
time_conv.push_back((tau[i]-0.5)/k2_inv*(h*h*1.0e-12)/Di[i]);
}
}
}
//read ion related list
//NOTE: ion diffusivity has INPUT unit: [m^2/sec]
// it must be converted to LB unit: [lu^2/lt]
if (ion_db->keyExists("IonDiffusivityList")){
IonDiffusivity.clear();
IonDiffusivity = ion_db->getVector<double>( "IonDiffusivityList" );
// time relaxation parameters tau also needs update
tau.clear();
if (IonDiffusivity.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
}
else{
for (int i=0; i<IonDiffusivity.size();i++){
//First, convert ion diffusivity in physical unit to LB unit
IonDiffusivity[i] = IonDiffusivity[i]*time_conv/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
//Second, re-calculate tau
tau.push_back(0.5+k2_inv*IonDiffusivity[i]);
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
}
}
else {
for (int i=0; i<IonDiffusivity.size();i++){
//convert ion diffusivity in physical unit to LB unit
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
}
// read time relaxation time list
//read ion algebric valence list
if (ion_db->keyExists("IonValenceList")){
IonValence.clear();
@@ -114,33 +148,255 @@ void ScaLBL_IonModel::ReadParams(string filename,int num_iter,int num_iter_Stoke
}
}
}
else {
for (int i=0; i<IonConcentration.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
//Read solid boundary condition specific to Ion model
BoundaryConditionSolid = 0;
if (ion_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = ion_db->getScalar<int>( "BC_Solid" );
}
// Read boundary condition for ion transport
// BC = 0: normal periodic BC
// BC = 1: fixed inlet and outlet ion concentration
BoundaryCondition = 0;
if (ion_db->keyExists( "BC" )){
BoundaryCondition = ion_db->getScalar<int>( "BC" );
}
if (BoundaryCondition==1){
//read boundary ion concentration list; INPUT unit [mol/m^3]
//it must be converted to LB unit [mol/lu^3]
//inlet
if (ion_db->keyExists("CinList")){
Cin.clear();
Cin = ion_db->getVector<double>( "CinList" );
if (Cin.size()!=number_ion_species){
ERROR("Error: number_ion_species and CinList must be the same length! \n");
}
else{
for (int i=0; i<Cin.size();i++){
Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
else {
for (int i=0; i<Cin.size();i++){
Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Ion Transport Solver: \n");
if (rank==0) printf(" Time conversion factor: %.5g [sec/lt]\n", time_conv);
if (rank==0) printf(" Internal iteration: %i [lt]\n", timestepMax);
for (int i=0; i<number_ion_species;i++){
if (rank==0) printf(" Ion %i: LB relaxation tau = %.5g\n", i+1,tau[i]);
}
//outlet
if (ion_db->keyExists("CoutList")){
Cout.clear();
Cout = ion_db->getVector<double>( "CoutList" );
if (Cout.size()!=number_ion_species){
ERROR("Error: number_ion_species and CoutList must be the same length! \n");
}
else{
for (int i=0; i<Cout.size();i++){
Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
else {
for (int i=0; i<Cout.size();i++){
Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
if (rank==0) printf("*****************************************************\n");
}
switch (BoundaryConditionSolid){
case 0:
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned\n");
break;
case 1:
if (rank==0) printf("LB Ion Solver: solid boundary: Neumann-type surfacen ion concentration is assigned\n");
break;
default:
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned\n");
break;
void ScaLBL_IonModel::ReadParams(string filename){
//NOTE: the maximum iteration timesteps for ions are left unspecified
// it relies on the multiphys controller to compute the max timestep
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
ion_db = db->getDatabase( "Ions" );
// Universal constant
kb = 1.38e-23;//Boltzmann constant;unit [J/K]
electron_charge = 1.6e-19;//electron charge;unit [C]
//---------------------- Default model parameters --------------------------//
T = 300.0;//temperature; unit [K]
Vt = kb*T/electron_charge;//thermal voltage; unit [Vy]
k2_inv = 4.0;//speed of sound for D3Q7 lattice
h = 1.0;//resolution; unit: um/lu
tolerance = 1.0e-8;
number_ion_species = 1;
tau.push_back(1.0);
IonDiffusivity.push_back(1.0e-9);//user-input diffusivity has physical unit [m^2/sec]
IonValence.push_back(1);//algebraic valence charge
IonConcentration.push_back(1.0e-3);//user-input ion concentration has physical unit [mol/m^3]
Cin.push_back(1.0e-3);//user-input inlet boundary ion concentration;unit [mol/m^3]
Cout.push_back(1.0e-3);//user-input outlet boundary ion concentration;unit [mol/m^3]
//tau.push_back(0.5+k2_inv*time_conv/(h*1.0e-6)/(h*1.0e-6)*IonDiffusivity[0]);
time_conv.push_back((tau[0]-0.5)/k2_inv*(h*h*1.0e-12)/IonDiffusivity[0]);
fluidVelx_dummy = 0.0;//for debugging, unit [m/sec]
fluidVely_dummy = 0.0;//for debugging, unit [m/sec]
fluidVelz_dummy = 0.0;//for debugging, unit [m/sec]
Ex_dummy = 0.0;//for debugging, unit [V/m]
Ey_dummy = 0.0;//for debugging, unit [V/m]
Ez_dummy = 0.0;//for debugging, unit [V/m]
//--------------------------------------------------------------------------//
// Read domain parameters
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
}
// LB-Ion Model parameters
//if (ion_db->keyExists( "timestepMax" )){
// timestepMax = ion_db->getScalar<int>( "timestepMax" );
//}
if (ion_db->keyExists( "tolerance" )){
tolerance = ion_db->getScalar<double>( "tolerance" );
}
if (ion_db->keyExists( "temperature" )){
T = ion_db->getScalar<int>( "temperature" );
//re-calculate thermal voltage
Vt = kb*T/electron_charge;//thermal voltage; unit [Vy]
}
if (ion_db->keyExists( "FluidVelDummy" )){
fluidVelx_dummy = ion_db->getVector<double>( "FluidVelDummy" )[0];
fluidVely_dummy = ion_db->getVector<double>( "FluidVelDummy" )[1];
fluidVelz_dummy = ion_db->getVector<double>( "FluidVelDummy" )[2];
}
if (ion_db->keyExists( "ElectricFieldDummy" )){
Ex_dummy = ion_db->getVector<double>( "ElectricFieldDummy" )[0];
Ey_dummy = ion_db->getVector<double>( "ElectricFieldDummy" )[1];
Ez_dummy = ion_db->getVector<double>( "ElectricFieldDummy" )[2];
}
if (ion_db->keyExists( "number_ion_species" )){
number_ion_species = ion_db->getScalar<int>( "number_ion_species" );
}
if (ion_db->keyExists("tauList")){
tau.clear();
tau = ion_db->getVector<double>( "tauList" );
vector<double>Di = ion_db->getVector<double>( "IonDiffusivityList" );//temp storing ion diffusivity in physical unit
if (tau.size()!=number_ion_species || Di.size()!=number_ion_species){
ERROR("Error: number_ion_species, tauList and IonDiffusivityList must be of the same length! \n");
}
else{
time_conv.clear();
for (int i=0; i<tau.size();i++){
time_conv.push_back((tau[i]-0.5)/k2_inv*(h*h*1.0e-12)/Di[i]);
}
}
}
//read ion related list
//NOTE: ion diffusivity has INPUT unit: [m^2/sec]
// it must be converted to LB unit: [lu^2/lt]
if (ion_db->keyExists("IonDiffusivityList")){
IonDiffusivity.clear();
IonDiffusivity = ion_db->getVector<double>( "IonDiffusivityList" );
if (IonDiffusivity.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n");
}
else{
for (int i=0; i<IonDiffusivity.size();i++){
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
}
}
else {
for (int i=0; i<IonDiffusivity.size();i++){
//convert ion diffusivity in physical unit to LB unit
IonDiffusivity[i] = IonDiffusivity[i]*time_conv[i]/(h*h*1.0e-12);//LB diffusivity has unit [lu^2/lt]
}
}
// read time relaxation time list
//read ion algebric valence list
if (ion_db->keyExists("IonValenceList")){
IonValence.clear();
IonValence = ion_db->getVector<int>( "IonValenceList" );
if (IonValence.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonValenceList must be the same length! \n");
}
}
//read initial ion concentration list; INPUT unit [mol/m^3]
//it must be converted to LB unit [mol/lu^3]
if (ion_db->keyExists("IonConcentrationList")){
IonConcentration.clear();
IonConcentration = ion_db->getVector<double>( "IonConcentrationList" );
if (IonConcentration.size()!=number_ion_species){
ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n");
}
else{
for (int i=0; i<IonConcentration.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
else {
for (int i=0; i<IonConcentration.size();i++){
IonConcentration[i] = IonConcentration[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
//Read solid boundary condition specific to Ion model
BoundaryConditionSolid = 0;
if (ion_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = ion_db->getScalar<int>( "BC_Solid" );
}
// Read boundary condition for ion transport
// BC = 0: normal periodic BC
// BC = 1: fixed inlet and outlet ion concentration
BoundaryCondition = 0;
if (ion_db->keyExists( "BC" )){
BoundaryCondition = ion_db->getScalar<int>( "BC" );
}
if (BoundaryCondition==1){
//read boundary ion concentration list; INPUT unit [mol/m^3]
//it must be converted to LB unit [mol/lu^3]
//inlet
if (ion_db->keyExists("CinList")){
Cin.clear();
Cin = ion_db->getVector<double>( "CinList" );
if (Cin.size()!=number_ion_species){
ERROR("Error: number_ion_species and CinList must be the same length! \n");
}
else{
for (int i=0; i<Cin.size();i++){
Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
else {
for (int i=0; i<Cin.size();i++){
Cin[i] = Cin[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
//outlet
if (ion_db->keyExists("CoutList")){
Cout.clear();
Cout = ion_db->getVector<double>( "CoutList" );
if (Cout.size()!=number_ion_species){
ERROR("Error: number_ion_species and CoutList must be the same length! \n");
}
else{
for (int i=0; i<Cout.size();i++){
Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
else {
for (int i=0; i<Cout.size();i++){
Cout[i] = Cout[i]*(h*h*h*1.0e-18);//LB ion concentration has unit [mol/lu^3]
}
}
}
}
@@ -261,7 +517,7 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
//NOTE need to convert the user input phys unit to LB unit
AFFINITY = AFFINITY*(h*h*1.0e-12);
AFFINITY = AFFINITY*(h*h*h*1.0e-18);
label_count[idx] += 1.0;
idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component
@@ -360,6 +616,27 @@ void ScaLBL_IonModel::Initialize(){
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
}
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Ion Transport Solver: \n");
for (int i=0; i<number_ion_species;i++){
if (rank==0) printf(" Ion %i: LB relaxation tau = %.5g\n", i+1,tau[i]);
if (rank==0) printf(" Time conversion factor: %.5g [sec/lt]\n", time_conv[i]);
if (rank==0) printf(" Internal iteration: %i [lt]\n", timestepMax[i]);
}
if (rank==0) printf("*****************************************************\n");
switch (BoundaryConditionSolid){
case 0:
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned\n");
break;
case 1:
if (rank==0) printf("LB Ion Solver: solid boundary: Dirichlet-type surfacen ion concentration is assigned\n");
break;
default:
if (rank==0) printf("LB Ion Solver: solid boundary: non-flux boundary is assigned\n");
break;
}
}
void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
@@ -379,77 +656,73 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//starttime = MPI_Wtime();
timestep=0;
while (timestep < timestepMax) {
//************************************************************************/
// *************ODD TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FROM NORMAL
for (int ic=0; ic<number_ion_species; ic++){
timestep=0;
while (timestep < timestepMax[ic]) {
//************************************************************************/
// *************ODD TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_IonConcentration(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
// Set boundary conditions
if (BoundaryCondition == 1){
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z(NeighborList, &fq[ic*Np*7], Cin[ic], timestep);
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z(NeighborList, &fq[ic*Np*7], Cout[ic], timestep);
}
//-------------------------//
ScaLBL_D3Q7_AAodd_IonConcentration(NeighborList, &fq[ic*Np*7],&Ci[ic*Np], 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
}
//LB-Ion collison
for (int ic=0; ic<number_ion_species; ic++){
//LB-Ion collison
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
// Set boundary conditions
/* ... */
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAodd_Ion(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
}
if (BoundaryConditionSolid==1){
for (int ic=0; ic<number_ion_species; ic++){
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidNeumannD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
}
}
// *************EVEN TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FORM NORMAL
if (BoundaryConditionSolid==1){
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
}
// *************EVEN TIMESTEP*************//
timestep++;
//Update ion concentration and charge density
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
// Set boundary conditions
if (BoundaryCondition == 1){
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_z(NeighborList, &fq[ic*Np*7], Cin[ic], timestep);
ScaLBL_Comm->D3Q7_Ion_Concentration_BC_Z(NeighborList, &fq[ic*Np*7], Cout[ic], timestep);
}
//-------------------------//
ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic*Np*7],&Ci[ic*Np], 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
}
//LB-Ion collison
for (int ic=0; ic<number_ion_species; ic++){
//LB-Ion collison
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
}
// Set boundary conditions
/* ... */
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_AAeven_Ion(&fq[ic*Np*7],&Ci[ic*Np],Velocity,ElectricField,IonDiffusivity[ic],IonValence[ic],
rlx[ic],Vt,0, ScaLBL_Comm->LastExterior(), Np);
}
if (BoundaryConditionSolid==1){
for (int ic=0; ic<number_ion_species; ic++){
if (BoundaryConditionSolid==1){
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidNeumannD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
}
}
//************************************************************************/
}
}
//Compute charge density for Poisson equation
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np);
}
//************************************************************************/
//stoptime = MPI_Wtime();
//if (rank==0) printf("-------------------------------------------------------------------\n");
@@ -464,21 +737,15 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//MLUPS *= nprocs;
//if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
//if (rank==0) printf("********************************************************\n");
}
//TODO this ruin the ion concentration on device
//need to do something similar to electric field
void ScaLBL_IonModel::getIonConcentration(int timestep){
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_IonConcentration_Phys(Ci, h, ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_IonConcentration_Phys(Ci, h, ic, 0, ScaLBL_Comm->LastExterior(), Np);
}
DoubleArray PhaseField(Nx,Ny,Nz);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
IonConcentration_LB_to_Phys(PhaseField);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
@@ -489,6 +756,100 @@ void ScaLBL_IonModel::getIonConcentration(int timestep){
}
void ScaLBL_IonModel::IonConcentration_LB_to_Phys(DoubleArray &Den_reg){
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int idx=Map(i,j,k);
if (!(idx < 0)){
Den_reg(i,j,k) = Den_reg(i,j,k)/(h*h*h*1.0e-18);
}
}
}
}
}
void ScaLBL_IonModel::DummyFluidVelocity(){
double *FluidVelocity_host;
FluidVelocity_host = new double[3*Np];
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
int idx=Map(i,j,k);
if (!(idx < 0))
FluidVelocity_host[idx+0*Np] = fluidVelx_dummy/(h*1.0e-6)*time_conv[0];
FluidVelocity_host[idx+1*Np] = fluidVely_dummy/(h*1.0e-6)*time_conv[0];
FluidVelocity_host[idx+2*Np] = fluidVelz_dummy/(h*1.0e-6)*time_conv[0];
}
}
}
ScaLBL_AllocateDeviceMemory((void **) &FluidVelocityDummy, sizeof(double)*3*Np);
ScaLBL_CopyToDevice(FluidVelocityDummy, FluidVelocity_host, sizeof(double)*3*Np);
ScaLBL_DeviceBarrier();
delete [] FluidVelocity_host;
}
void ScaLBL_IonModel::DummyElectricField(){
double *ElectricField_host;
ElectricField_host = new double[3*Np];
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
int idx=Map(i,j,k);
if (!(idx < 0))
ElectricField_host[idx+0*Np] = Ex_dummy*(h*1.0e-6);
ElectricField_host[idx+1*Np] = Ey_dummy*(h*1.0e-6);
ElectricField_host[idx+2*Np] = Ez_dummy*(h*1.0e-6);
}
}
}
ScaLBL_AllocateDeviceMemory((void **) &ElectricFieldDummy, sizeof(double)*3*Np);
ScaLBL_CopyToDevice(ElectricFieldDummy, ElectricField_host, sizeof(double)*3*Np);
ScaLBL_DeviceBarrier();
delete [] ElectricField_host;
}
double ScaLBL_IonModel::CalIonDenConvergence(vector<double> &ci_avg_previous){
double *Ci_host;
Ci_host = new double[Np];
vector<double> error(number_ion_species,0.0);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_CopyToHost(Ci_host,&Ci[ic*Np],Np*sizeof(double));
double count_loc=0;
double count;
double ci_avg;
double ci_loc=0.f;
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
ci_loc +=Ci_host[idx];
count_loc+=1.0;
}
for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
ci_loc +=Ci_host[idx];
count_loc+=1.0;
}
MPI_Allreduce(&ci_loc,&ci_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
ci_avg /= count;
double ci_avg_mag=ci_avg;
if (ci_avg==0.0) ci_avg_mag=1.0;
error[ic] = fabs(ci_avg-ci_avg_previous[ic])/fabs(ci_avg_mag);
ci_avg_previous[ic] = ci_avg;
}
double error_max;
error_max = *max_element(error.begin(),error.end());
if (rank==0){
printf("IonModel: error max: %.5g\n",error_max);
}
return error_max;
}
//void ScaLBL_IonModel::getIonConcentration(){
// for (int ic=0; ic<number_ion_species; ic++){
// ScaLBL_IonConcentration_Phys(Ci, h, ic, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);

20
models/IonModel.h
View File

@@ -22,7 +22,8 @@ public:
~ScaLBL_IonModel();
// functions in they should be run
void ReadParams(string filename,int num_iter,int num_iter_Stokes,double time_conv_Stokes);
void ReadParams(string filename,vector<int> &num_iter);
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void ReadInput();
@@ -30,24 +31,30 @@ public:
void Initialize();
void Run(double *Velocity, double *ElectricField);
void getIonConcentration(int timestep);
void DummyFluidVelocity();
void DummyElectricField();
double CalIonDenConvergence(vector<double> &ci_avg_previous);
//bool Restart,pBC;
int timestep,timestepMax;
int timestep;
vector<int> timestepMax;
int BoundaryCondition;
int BoundaryConditionSolid;
double h;//domain resolution, unit [um/lu]
double time_conv;
double kb,electron_charge,T,Vt;
double k2_inv;
double tolerance;
double Ex,Ey,Ez;
double fluidVelx_dummy,fluidVely_dummy,fluidVelz_dummy;
double Ex_dummy,Ey_dummy,Ez_dummy;
int number_ion_species;
vector<double> IonDiffusivity;//User input unit [m^2/sec]
vector<int> IonValence;
vector<double> IonConcentration;//unit [mol/m^3]
//vector<double> deltaT;
vector<double> Cin;//unit [mol/m^3]
vector<double> Cout;//unit [mol/m^3]
vector<double> tau;
vector<double> time_conv;
int Nx,Ny,Nz,N,Np;
int rank,nprocx,nprocy,nprocz,nprocs;
@@ -68,6 +75,8 @@ public:
double *Ci;
double *ChargeDensity;
double *IonSolid;
double *FluidVelocityDummy;
double *ElectricFieldDummy;
private:
MPI_Comm comm;
@@ -80,4 +89,5 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignSolidBoundary(double *ion_solid);
void IonConcentration_LB_to_Phys(DoubleArray &Den_reg);
};

119
models/MultiPhysController.cpp
View File

@@ -1,7 +1,8 @@
#include "models/MultiPhysController.h"
ScaLBL_Multiphys_Controller::ScaLBL_Multiphys_Controller(int RANK, int NP, MPI_Comm COMM):
rank(RANK),nprocs(NP),Restart(0),timestepMax(0),num_iter_Stokes(0),num_iter_Ion(0),SchmidtNum(0),comm(COMM)
rank(RANK),nprocs(NP),Restart(0),timestepMax(0),num_iter_Stokes(0),num_iter_Ion(0),
analysis_interval(0),tolerance(0),comm(COMM)
{
}
@@ -19,34 +20,50 @@ void ScaLBL_Multiphys_Controller::ReadParams(string filename){
// Default parameters
timestepMax = 10000;
Restart = false;
SchmidtNum = 1.0;
num_iter_Stokes=1;
num_iter_Ion=1;
num_iter_Ion.push_back(1);
analysis_interval = 500;
tolerance = 1.0e-6;
// load input parameters
if (study_db->keyExists( "timestepMax" )){
timestepMax = study_db->getScalar<int>( "timestepMax" );
}
if (study_db->keyExists( "Schmidt_Number" )){
SchmidtNum = study_db->getScalar<double>( "Schmidt_Number" );
if (study_db->keyExists( "analysis_interval" )){
analysis_interval = study_db->getScalar<int>( "analysis_interval" );
}
if (study_db->keyExists( "tolerance" )){
tolerance = study_db->getScalar<double>( "tolerance" );
}
//if (study_db->keyExists( "time_conv" )){
// time_conv = study_db->getScalar<double>( "time_conv" );
//}
//if (study_db->keyExists( "Schmidt_Number" )){
// SchmidtNum = study_db->getScalar<double>( "Schmidt_Number" );
//}
// recalculate relevant parameters
if (SchmidtNum>1){
num_iter_Stokes = int(round(SchmidtNum/2)*2);
num_iter_Ion = 1;
}
else if (SchmidtNum>0 && SchmidtNum<1){
num_iter_Ion = int(round((1.0/SchmidtNum)/2)*2);
num_iter_Stokes = 1;
}
else{
ERROR("Error: SchmidtNum (Schmidt number) must be a positive number! \n");
}
//if (SchmidtNum>1){
// num_iter_Stokes = int(round(SchmidtNum/2)*2);
// num_iter_Ion = 1;
//}
//else if (SchmidtNum>0 && SchmidtNum<1){
// num_iter_Ion = int(round((1.0/SchmidtNum)/2)*2);
// num_iter_Stokes = 1;
//}
//else if (SchmidtNum==1){
// num_iter_Stokes = 1;
// num_iter_Ion = 1;
//}
//else{
// ERROR("Error: SchmidtNum (Schmidt number) must be a positive number! \n");
//}
// load input parameters
// in case user wants to have an absolute control over the iternal iteration
if (study_db->keyExists( "num_iter_Ion" )){
num_iter_Ion = study_db->getScalar<int>( "num_iter_Ion" );
if (study_db->keyExists( "num_iter_Ion_List" )){
num_iter_Ion.clear();
num_iter_Ion = study_db->getVector<int>( "num_iter_Ion_List" );
}
if (study_db->keyExists( "num_iter_Stokes" )){
num_iter_Stokes = study_db->getScalar<int>( "num_iter_Stokes" );
@@ -54,4 +71,70 @@ void ScaLBL_Multiphys_Controller::ReadParams(string filename){
}
int ScaLBL_Multiphys_Controller::getStokesNumIter_PNP_coupling(double StokesTimeConv,const vector<double> &IonTimeConv){
//Return number of internal iterations for the Stokes solver
int num_iter_stokes;
vector<double> TimeConv;
printf("*****Debug; IonTimeConv size = %i\n",IonTimeConv.size());
for (unsigned int i =0; i<IonTimeConv.size();i++){
printf("*****Debug; Ion %i; IonTimeConv = %.5g\n",i,IonTimeConv[i]);
}
TimeConv.assign(IonTimeConv.begin(),IonTimeConv.end());
TimeConv.insert(TimeConv.begin(),StokesTimeConv);
for (unsigned int i =0; i<TimeConv.size();i++){
printf("*****Debug; all TimeConv %i; TimeConv = %.5g\n",i,TimeConv[i]);
}
vector<double>::iterator it_max = max_element(TimeConv.begin(),TimeConv.end());
int idx_max = distance(TimeConv.begin(),it_max);
if (idx_max==0){
num_iter_stokes = 2;
}
else{
double temp = 2*TimeConv[idx_max]/StokesTimeConv;//the factor 2 is the number of iterations for the element has max time_conv
num_iter_stokes = int(round(temp/2)*2);
}
return num_iter_stokes;
}
vector<int> ScaLBL_Multiphys_Controller::getIonNumIter_PNP_coupling(double StokesTimeConv,const vector<double> &IonTimeConv){
//Return number of internal iterations for the Ion transport solver
vector<int> num_iter_ion;
vector<double> TimeConv;
TimeConv.assign(IonTimeConv.begin(),IonTimeConv.end());
TimeConv.insert(TimeConv.begin(),StokesTimeConv);
vector<double>::iterator it_max = max_element(TimeConv.begin(),TimeConv.end());
unsigned int idx_max = distance(TimeConv.begin(),it_max);
if (idx_max==0){
for (unsigned int idx=1;idx<TimeConv.size();idx++){
double temp = 2*StokesTimeConv/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==1){
num_iter_ion.push_back(2);
for (unsigned int idx=2;idx<TimeConv.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==TimeConv.size()-1){
for (unsigned int idx=1;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=1;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<TimeConv.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;
}

12
models/MultiPhysController.h
View File

@@ -8,6 +8,8 @@
#include <exception>
#include <stdexcept>
#include <fstream>
#include <vector>
#include <algorithm>
#include "common/ScaLBL.h"
#include "common/Communication.h"
@@ -22,13 +24,17 @@ public:
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
int getStokesNumIter_PNP_coupling(double StokesTimeConv,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);
bool Restart;
//int timestep;
int timestepMax;
int num_iter_Stokes;
int num_iter_Ion;
double SchmidtNum;//Schmidt number = kinematic_viscosity/mass_diffusivity
vector<int> num_iter_Ion;
int analysis_interval;
double tolerance;
//double SchmidtNum;//Schmidt number = kinematic_viscosity/mass_diffusivity
int rank,nprocs;

2
models/PoissonSolver.cpp
View File

@@ -59,7 +59,7 @@ void ScaLBL_Poisson::ReadParams(string filename){
if (electric_db->keyExists( "BC_Solid" )){
BoundaryConditionSolid = electric_db->getScalar<int>( "BC_Solid" );
}
// Read boundary condition for electric potentiona
// Read boundary condition for electric potential
// BC = 0: normal periodic BC
// BC = 1: fixed inlet and outlet potential
BoundaryCondition = 0;

234
models/StokesModel.cpp
View File

@@ -91,12 +91,84 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){
time_conv = (h*h*1.0e-12)*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: \n");
if (rank==0) printf(" Time conversion factor: %.5g [sec/lt]\n", time_conv);
if (rank==0) printf(" Internal iteration: %i [lt]\n", timestepMax);
if (rank==0) printf("*****************************************************\n");
}
void ScaLBL_StokesModel::ReadParams(string filename){
//NOTE the max time step is left unspecified
// read the input database
db = std::make_shared<Database>( filename );
domain_db = db->getDatabase( "Domain" );
stokes_db = db->getDatabase( "Stokes" );
//---------------------- Default model parameters --------------------------//
rho_phys = 1000.0; //by default use water density; unit [kg/m^3]
nu_phys = 1.004e-6;//by default use water kinematic viscosity at 20C; unit [m^2/sec]
h = 1.0;//image resolution;[um]
tau = 1.0;
mu = (tau-0.5)/3.0;//LB kinematic viscosity;unit [lu^2/lt]
time_conv = h*h*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
rho0 = 1.0;//LB density
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
//--------------------------------------------------------------------------//
// Read domain parameters
BoundaryCondition = 0;
if (domain_db->keyExists( "BC" )){
BoundaryCondition = domain_db->getScalar<int>( "BC" );
}
if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu
h = domain_db->getScalar<double>( "voxel_length" );
}
// Single-fluid Navier-Stokes Model parameters
//if (stokes_db->keyExists( "timestepMax" )){
// timestepMax = stokes_db->getScalar<int>( "timestepMax" );
//}
if (stokes_db->keyExists( "tolerance" )){
tolerance = stokes_db->getScalar<double>( "tolerance" );
}
if (stokes_db->keyExists( "tau" )){
tau = stokes_db->getScalar<double>( "tau" );
}
if (stokes_db->keyExists( "rho0" )){
rho0 = stokes_db->getScalar<double>( "rho0" );
}
if (stokes_db->keyExists( "nu_phys" )){
nu_phys = stokes_db->getScalar<double>( "nu_phys" );
}
if (stokes_db->keyExists( "rho_phys" )){
rho_phys = stokes_db->getScalar<double>( "rho_phys" );
}
if (stokes_db->keyExists( "F" )){
Fx = stokes_db->getVector<double>( "F" )[0];
Fy = stokes_db->getVector<double>( "F" )[1];
Fz = stokes_db->getVector<double>( "F" )[2];
}
if (stokes_db->keyExists( "Restart" )){
Restart = stokes_db->getScalar<bool>( "Restart" );
}
if (stokes_db->keyExists( "din" )){
din = stokes_db->getScalar<double>( "din" );
}
if (stokes_db->keyExists( "dout" )){
dout = stokes_db->getScalar<double>( "dout" );
}
if (stokes_db->keyExists( "flux" )){
flux = stokes_db->getScalar<double>( "flux" );
}
// Re-calculate model parameters due to parameter read
mu=(tau-0.5)/3.0;
time_conv = (h*h*1.0e-12)*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt]
den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density
}
void ScaLBL_StokesModel::SetDomain(){
Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis
Mask = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // mask domain removes immobile phases
@@ -235,6 +307,12 @@ void ScaLBL_StokesModel::Initialize(){
if (rank==0) printf("LB Single-Fluid Solver: Initializing distributions \n");
if (rank==0) printf("****************************************************************\n");
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf("*****************************************************\n");
if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: \n");
if (rank==0) printf(" Time conversion factor: %.5g [sec/lt]\n", time_conv);
if (rank==0) printf(" Internal iteration: %i [lt]\n", timestepMax);
if (rank==0) printf("*****************************************************\n");
}
void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
@@ -243,6 +321,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
timestep = 0;
while (timestep < timestepMax) {
//************************************************************************/
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@@ -334,34 +413,123 @@ void ScaLBL_StokesModel::Velocity_LB_to_Phys(DoubleArray &Vel_reg){
}
}
//void ScaLBL_StokesModel::getVelocity(){
// //get velocity in physical unit [m/sec]
// ScaLBL_D3Q19_Momentum_Phys(fq, Velocity, h, time_conv, Np);
// ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//
// DoubleArray PhaseField(Nx,Ny,Nz);
// ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
// FILE *VELX_FILE;
// sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
// VELX_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,VELX_FILE);
// fclose(VELX_FILE);
//
// ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
// FILE *VELY_FILE;
// sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
// VELY_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,VELY_FILE);
// fclose(VELY_FILE);
//
// ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
// FILE *VELZ_FILE;
// sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
// VELZ_FILE = fopen(LocalRankFilename,"wb");
// fwrite(PhaseField.data(),8,N,VELZ_FILE);
// fclose(VELZ_FILE);
//
//}
vector<double> ScaLBL_StokesModel::computeElectricForceAvg(double *ChargeDensity, double *ElectricField){
double *Ex_host;
double *Ey_host;
double *Ez_host;
Ex_host = new double[Np];
Ey_host = new double[Np];
Ez_host = new double[Np];
double *rhoE_host;
rhoE_host = new double[Np];
ScaLBL_CopyToHost(Ex_host,&ElectricField[0*Np],Np*sizeof(double));
ScaLBL_CopyToHost(Ey_host,&ElectricField[1*Np],Np*sizeof(double));
ScaLBL_CopyToHost(Ez_host,&ElectricField[2*Np],Np*sizeof(double));
ScaLBL_CopyToHost(rhoE_host,ChargeDensity,Np*sizeof(double));
double count_loc=0;
double count;
double Fx_avg,Fy_avg,Fz_avg;//average electric field induced force
double Fx_loc,Fy_loc,Fz_loc;
Fx_loc = Fy_loc = Fz_loc = 0.0;
for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
Fx_loc += rhoE_host[idx]*Ex_host[idx]*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fy_loc += rhoE_host[idx]*Ey_host[idx]*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz_loc += rhoE_host[idx]*Ez_host[idx]*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
count_loc+=1.0;
}
for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
Fx_loc += rhoE_host[idx]*Ex_host[idx]*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fy_loc += rhoE_host[idx]*Ey_host[idx]*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
Fz_loc += rhoE_host[idx]*Ez_host[idx]*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
count_loc+=1.0;
}
MPI_Allreduce(&Fx_loc,&Fx_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&Fy_loc,&Fy_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&Fz_loc,&Fz_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
Fx_avg /= count;
Fy_avg /= count;
Fz_avg /= count;
vector<double>F_avg{Fx_avg,Fy_avg,Fz_avg};
delete [] Ex_host;
delete [] Ey_host;
delete [] Ez_host;
delete [] rhoE_host;
return F_avg;
}
double ScaLBL_StokesModel::CalVelocityConvergence(double& flow_rate_previous,double *ChargeDensity, double *ElectricField){
//-----------------------------------------------------
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
double count_loc=0;
double count;
double vax,vay,vaz;
double vax_loc,vay_loc,vaz_loc;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
if (Distance(i,j,k) > 0){
vax_loc += Velocity_x(i,j,k);
vay_loc += Velocity_y(i,j,k);
vaz_loc += Velocity_z(i,j,k);
count_loc+=1.0;
}
}
}
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax /= count;
vay /= count;
vaz /= count;
vector<double> Eforce;
Eforce = computeElectricForceAvg(ChargeDensity,ElectricField);
double TFx = Fx+Eforce[0];//TF: total body force
double TFy = Fy+Eforce[1];
double TFz = Fz+Eforce[2];
double force_mag = sqrt(TFx*TFx+TFy*TFy+TFz*TFz);
double dir_x = TFx/force_mag;
double dir_y = TFy/force_mag;
double dir_z = TFz/force_mag;
if (force_mag == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
force_mag = 1.0;
}
double flow_rate = (vax*dir_x + vay*dir_y + vaz*dir_z);
double error = fabs(flow_rate - flow_rate_previous) / fabs(flow_rate);
flow_rate_previous = flow_rate;
//----------------------------------------------------
//for debugging
if (rank==0){
printf("StokesModel: error: %.5g\n",error);
}
return error;
}
void ScaLBL_StokesModel::Run(){
double rlx_setA=1.0/tau;

3
models/StokesModel.h
View File

@@ -22,6 +22,7 @@ public:
// functions in they should be run
void ReadParams(string filename,int num_iter);
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void SetDomain();
void ReadInput();
@@ -31,6 +32,7 @@ public:
void Run_Lite(double *ChargeDensity, double *ElectricField);
void VelocityField();
void getVelocity(int timestep);
double CalVelocityConvergence(double& flow_rate_previous,double *ChargeDensity, double *ElectricField);
bool Restart,pBC;
int timestep,timestepMax;
@@ -81,4 +83,5 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void Velocity_LB_to_Phys(DoubleArray &Vel_reg);
vector<double> computeElectricForceAvg(double *ChargeDensity, double *ElectricField);
};

3
tests/CMakeLists.txt
View File

@@ -48,6 +48,9 @@ ADD_LBPM_TEST( TestTopo3D )
ADD_LBPM_TEST( TestFluxBC )
ADD_LBPM_TEST( TestMap )
ADD_LBPM_TEST( TestPoissonSolver )
ADD_LBPM_TEST( TestIonModel )
ADD_LBPM_TEST( TestNernstPlanck )
ADD_LBPM_TEST( TestPNP_Stokes )
#ADD_LBPM_TEST( TestMRT )
#ADD_LBPM_TEST( TestColorGrad )
#ADD_LBPM_TEST( TestColorGradDFH )

89
tests/TestIonModel.cpp Normal file
View File

@@ -0,0 +1,89 @@
#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/MultiPhysController.h"
using namespace std;
//***************************************************************************
// Test lattice-Boltzmann Ion Model coupled with Poisson equation
//***************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI
int provided_thread_support = -1;
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
MPI_Comm comm;
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
int rank = comm_rank(comm);
int nprocs = comm_size(comm);
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
}
// Limit scope so variables that contain communicators will free before MPI_Finialize
{
if (rank == 0){
printf("**************************************\n");
printf("Running Test for Ion Transport \n");
printf("**************************************\n");
}
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
auto filename = argv[1];
ScaLBL_IonModel IonModel(rank,nprocs,comm);
ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling
// Load controller information
Study.ReadParams(filename);
// Initialize LB-Ion model
IonModel.ReadParams(filename,Study.num_iter_Ion);
IonModel.SetDomain();
IonModel.ReadInput();
IonModel.Create();
IonModel.Initialize();
IonModel.DummyFluidVelocity();
IonModel.DummyElectricField();
int timestep=0;
double error = 1.0;
vector<double>ci_avg_previous{0.0,0.0};//assuming 1:1 solution
while (timestep < Study.timestepMax && error > Study.tolerance){
timestep++;
IonModel.Run(IonModel.FluidVelocityDummy,IonModel.ElectricFieldDummy); //solve for ion transport and electric potential
timestep++;//AA operations
if (timestep%Study.analysis_interval==0){
error = IonModel.CalIonDenConvergence(ci_avg_previous);
}
}
IonModel.getIonConcentration(timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("*************************************************************\n");
PROFILE_STOP("Main");
PROFILE_SAVE("TestIonModel",1);
// ****************************************************
MPI_Barrier(comm);
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Comm_free(&comm);
MPI_Finalize();
}

101
tests/TestNernstPlanck.cpp Normal file
View File

@@ -0,0 +1,101 @@
#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"
using namespace std;
//***************************************************************************
// Test lattice-Boltzmann Ion Model coupled with Poisson equation
//***************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI
int provided_thread_support = -1;
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
MPI_Comm comm;
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
int rank = comm_rank(comm);
int nprocs = comm_size(comm);
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
}
// Limit scope so variables that contain communicators will free before MPI_Finialize
{
if (rank == 0){
printf("********************************************************\n");
printf("Running Test for LB-Poisson-Ion Coupling \n");
printf("********************************************************\n");
}
//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);
// Initialize LB-Ion model
IonModel.ReadParams(filename,Study.num_iter_Ion);
IonModel.SetDomain();
IonModel.ReadInput();
IonModel.Create();
IonModel.Initialize();
IonModel.DummyFluidVelocity();
// Initialize LB-Poisson model
PoissonSolver.ReadParams(filename);
PoissonSolver.SetDomain();
PoissonSolver.ReadInput();
PoissonSolver.Create();
PoissonSolver.Initialize();
int timestep=0;
double error = 1.0;
vector<double>ci_avg_previous{0.0,0.0};//assuming 1:1 solution
while (timestep < Study.timestepMax && error > Study.tolerance){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity);//solve Poisson equtaion to get steady-state electrical potental
IonModel.Run(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField); //solve for ion transport and electric potential
timestep++;//AA operations
if (timestep%Study.analysis_interval==0){
error = IonModel.CalIonDenConvergence(ci_avg_previous);
}
}
PoissonSolver.getElectricPotential(timestep);
PoissonSolver.getElectricField(timestep);
IonModel.getIonConcentration(timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("*************************************************************\n");
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_electrokinetic_simulator",1);
// ****************************************************
MPI_Barrier(comm);
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Comm_free(&comm);
MPI_Finalize();
}

123
tests/TestPNP_Stokes.cpp Normal file
View File

@@ -0,0 +1,123 @@
#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/StokesModel.h"
#include "models/PoissonSolver.h"
#include "models/MultiPhysController.h"
using namespace std;
//***************************************************************************
// Test lattice-Boltzmann Ion Model coupled with Poisson equation
//***************************************************************************
int main(int argc, char **argv)
{
// Initialize MPI
int provided_thread_support = -1;
MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
MPI_Comm comm;
MPI_Comm_dup(MPI_COMM_WORLD,&comm);
int rank = comm_rank(comm);
int nprocs = comm_size(comm);
if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
}
// Limit scope so variables that contain communicators will free before MPI_Finialize
{
if (rank == 0){
printf("********************************************************\n");
printf("Running Test for LB-Poisson-Ion Coupling \n");
printf("********************************************************\n");
}
//PROFILE_ENABLE_TRACE();
//PROFILE_ENABLE_MEMORY();
PROFILE_SYNCHRONIZE();
PROFILE_START("Main");
Utilities::setErrorHandlers();
auto filename = argv[1];
ScaLBL_StokesModel StokesModel(rank,nprocs,comm);
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 Navier-Stokes and Ion solvers
StokesModel.ReadParams(filename);
IonModel.ReadParams(filename);
// Setup other model specific structures
StokesModel.SetDomain();
StokesModel.ReadInput();
StokesModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
IonModel.SetDomain();
IonModel.ReadInput();
IonModel.Create();
// Get internal iteration number
StokesModel.timestepMax = Study.getStokesNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
StokesModel.Initialize(); // initializing the model will set initial conditions for variables
IonModel.timestepMax = Study.getIonNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv);
IonModel.Initialize();
// Initialize LB-Poisson model
PoissonSolver.ReadParams(filename);
PoissonSolver.SetDomain();
PoissonSolver.ReadInput();
PoissonSolver.Create();
PoissonSolver.Initialize();
int timestep=0;
double error = 1.0;
double error_ion = 1.0;
double error_stokes = 1.0;
vector<double>ci_avg_previous{0.0,0.0};//assuming 1:1 solution
double vel_avg_previous = 0.0;
while (timestep < Study.timestepMax && error > Study.tolerance){
timestep++;
PoissonSolver.Run(IonModel.ChargeDensity);//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(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
timestep++;//AA operations
if (timestep%Study.analysis_interval==0){
error_ion = IonModel.CalIonDenConvergence(ci_avg_previous);
error_stokes = StokesModel.CalVelocityConvergence(vel_avg_previous,IonModel.ChargeDensity,PoissonSolver.ElectricField);
error = max(error_ion,error_stokes);
}
}
PoissonSolver.getElectricPotential(timestep);
PoissonSolver.getElectricField(timestep);
IonModel.getIonConcentration(timestep);
StokesModel.getVelocity(timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("*************************************************************\n");
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_electrokinetic_simulator",1);
// ****************************************************
MPI_Barrier(comm);
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Comm_free(&comm);
MPI_Finalize();
}

29
tests/TestPoissonSolver.cpp
View File

@@ -50,7 +50,6 @@ int main(int argc, char **argv)
PoissonSolver.ReadInput();
PoissonSolver.Create();
PoissonSolver.Initialize();
PoissonSolver.getElectricPotential(0);
//Initialize dummy charge density for test
PoissonSolver.DummyChargeDensity();
@@ -59,34 +58,6 @@ int main(int argc, char **argv)
PoissonSolver.getElectricPotential(1);
PoissonSolver.getElectricField(1);
//int timestep=0;
//while (timestep < Study.timestepMax){
//
// timestep++;
// //if (rank==0) printf("timestep=%i; running Poisson solver\n",timestep);
// PoissonSolver.Run(IonModel.ChargeDensity);//solve Poisson equtaion to get steady-state electrical potental
// //PoissonSolver.getElectricPotential(timestep);
// //if (rank==0) printf("timestep=%i; running StokesModel\n",timestep);
// StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
// //StokesModel.getVelocity(timestep);
// //if (rank==0) printf("timestep=%i; running Ion model\n",timestep);
// IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
// //IonModel.getIonConcentration(timestep);
//
//
// timestep++;//AA operations
// //--------------------------------------------
// //potentially leave analysis module for future
// //--------------------------------------------
//}
//StokesModel.getVelocity(timestep);
//PoissonSolver.getElectricPotential(timestep);
//PoissonSolver.getElectricField(timestep);
//IonModel.getIonConcentration(timestep);
if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
if (rank==0) printf("*************************************************************\n");

2
tests/lbpm_electrokinetic_dfh_simulator.cpp
View File

@@ -61,7 +61,7 @@ int main(int argc, char **argv)
StokesModel.Initialize(); // initializing the model will set initial conditions for variables
// Initialize LB-Ion model
IonModel.ReadParams(filename,Study.num_iter_Ion,Study.num_iter_Stokes,StokesModel.time_conv);
IonModel.ReadParams(filename,Study.num_iter_Ion);
IonModel.SetDomain();
IonModel.ReadInput();
IonModel.Create();