debug charge density problems on crusher

This commit is contained in:
James E McClure 2022-04-24 15:55:23 -04:00
parent 9043751281
commit 3642a6ae9b
5 changed files with 38 additions and 11 deletions

View File

@ -301,7 +301,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Ion(double *dist, double *Den, double *FluxDi
extern "C" void ScaLBL_D3Q7_Ion_Init(double *dist, double *Den, double DenInit, int Np);
extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, 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_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np);
// LBM Poisson solver

View File

@ -566,7 +566,7 @@ __global__ void dvc_ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, in
}
}
__global__ void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
int n;
double Ci;//ion concentration of species i
@ -579,10 +579,17 @@ __global__ void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDe
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if (n<finish) {
Ci = Den[n+ion_component*Np];
CD = ChargeDensity[n];
if (ion_component == 0) CD=0.0;
CD_tmp = F*IonValence*Ci;
ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
ChargeDensity[n] = CD + CD_tmp;
// Ci = Den[n+ion_component*Np];
// CD = ChargeDensity[n];
// CD_tmp = F*IonValence*Ci;
// ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
}
}
}
@ -660,7 +667,7 @@ extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np)
//cudaProfilerStop();
}
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_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_Ion_ChargeDensity<<<NBLOCKS,NTHREADS >>>(Den,ChargeDensity,IonValence,ion_component,start,finish,Np);

View File

@ -564,7 +564,7 @@ __global__ void dvc_ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, in
}
}
__global__ void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, int IonValence, int ion_component, int start, int finish, int Np){
__global__ void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
int n;
double Ci;//ion concentration of species i
@ -579,8 +579,9 @@ __global__ void dvc_ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDe
if (n<finish) {
Ci = Den[n+ion_component*Np];
CD = ChargeDensity[n];
if (ion_component == 0) CD=0.0;
CD_tmp = F*IonValence*Ci;
ChargeDensity[n] = CD*(ion_component>0) + CD_tmp;
ChargeDensity[n] = CD + CD_tmp;
}
}
}
@ -658,7 +659,7 @@ extern "C" void ScaLBL_D3Q7_Ion_Init_FromFile(double *dist, double *Den, int Np)
//cudaProfilerStop();
}
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_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity, double IonValence, int ion_component, int start, int finish, int Np){
//cudaProfilerStart();
dvc_ScaLBL_D3Q7_Ion_ChargeDensity<<<NBLOCKS,NTHREADS >>>(Den,ChargeDensity,IonValence,ion_component,start,finish,Np);

View File

@ -1443,17 +1443,35 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
//Compute charge density for Poisson equation
for (size_t ic = 0; ic < number_ion_species; ic++) {
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic,
int Valence = IonValence[ic];
if (rank==0) printf("compute charge density for ion %i, Valence =%i \n", ic,Valence);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic,
ScaLBL_Comm->FirstInterior(),
ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0,
ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, Valence, ic, 0,
ScaLBL_Comm->LastExterior(), Np);
}
DoubleArray Charge(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map, ChargeDensity, Charge);
double charge_sum=0.0;
double charge_sum_total=0.0;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
charge_sum += Charge(i,j,k);
}
}
}
printf(" Local charge value = %.8g (rank=%i)\n",charge_sum, rank);
ScaLBL_Comm->Barrier();
comm.barrier();
if (rank==0) printf(" IonMembrane: completeted full step \n");
fflush(stdout);
ScaLBL_Comm->Barrier();
comm.barrier();
//if (rank==0) printf(" IonMembrane: completeted full step \n");
//fflush(stdout);
//************************************************************************/
//if (rank==0) printf("-------------------------------------------------------------------\n");
//// Compute the walltime per timestep

View File

@ -557,6 +557,7 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
// Check convergence of steady-state solution
if (timestep==2){
//save electric potential for convergence check