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

initialize a new convergence scheme for Poisson solver

Browse Source
This commit is contained in:
Rex Zhe Li 2021-09-09 15:48:25 +10:00
parent a2ed3fdb34
commit c74b803157
4 changed files with 102 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
common/ScaLBL.h
View File

@ -132,6 +132,8 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q7_PoissonResidualError(int *neighborList, int *Map, double *ResidualError, double *Psi, double *Den_charge, double epsilon_LB,int strideY, int strideZ,int start, int finish);
//maybe deprecated
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
// int strideY, int strideZ,int start, int finish, int Np);

81
cpu/Poisson.cpp
View File

@ -235,6 +235,87 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in
}
}
extern "C" void ScaLBL_D3Q7_PoissonResidualError(int *neighborList, int *Map, double *ResidualError, double *Psi, double *Den_charge, double epsilon_LB,int strideY, int strideZ,int start, int finish){
int n,nn,ijk;
double psi;//electric potential
double rho_e;//local charge density
// neighbors of electric potential psi
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
double m0,m3,m5,m7;
double psi_Laplacian;
double residual_error;
for (n=start; n<finish; n++){
//Load data
rho_e = Den_charge[n];
ijk=Map[n];
psi = Psi[ijk];
// COMPUTE THE COLOR GRADIENT
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk-1; // neighbor index (get convention)
m1 = Psi[nn]; // get neighbor for phi - 1
//........................................................................
nn = ijk+1; // neighbor index (get convention)
m2 = Psi[nn]; // get neighbor for phi - 2
//........................................................................
nn = ijk-strideY; // neighbor index (get convention)
m3 = Psi[nn]; // get neighbor for phi - 3
//........................................................................
nn = ijk+strideY; // neighbor index (get convention)
m4 = Psi[nn]; // get neighbor for phi - 4
//........................................................................
nn = ijk-strideZ; // neighbor index (get convention)
m5 = Psi[nn]; // get neighbor for phi - 5
//........................................................................
nn = ijk+strideZ; // neighbor index (get convention)
m6 = Psi[nn]; // get neighbor for phi - 6
//........................................................................
nn = ijk-strideY-1; // neighbor index (get convention)
m7 = Psi[nn]; // get neighbor for phi - 7
//........................................................................
nn = ijk+strideY+1; // neighbor index (get convention)
m8 = Psi[nn]; // get neighbor for phi - 8
//........................................................................
nn = ijk+strideY-1; // neighbor index (get convention)
m9 = Psi[nn]; // get neighbor for phi - 9
//........................................................................
nn = ijk-strideY+1; // neighbor index (get convention)
m10 = Psi[nn]; // get neighbor for phi - 10
//........................................................................
nn = ijk-strideZ-1; // neighbor index (get convention)
m11 = Psi[nn]; // get neighbor for phi - 11
//........................................................................
nn = ijk+strideZ+1; // neighbor index (get convention)
m12 = Psi[nn]; // get neighbor for phi - 12
//........................................................................
nn = ijk+strideZ-1; // neighbor index (get convention)
m13 = Psi[nn]; // get neighbor for phi - 13
//........................................................................
nn = ijk-strideZ+1; // neighbor index (get convention)
m14 = Psi[nn]; // get neighbor for phi - 14
//........................................................................
nn = ijk-strideZ-strideY; // neighbor index (get convention)
m15 = Psi[nn]; // get neighbor for phi - 15
//........................................................................
nn = ijk+strideZ+strideY; // neighbor index (get convention)
m16 = Psi[nn]; // get neighbor for phi - 16
//........................................................................
nn = ijk+strideZ-strideY; // neighbor index (get convention)
m17 = Psi[nn]; // get neighbor for phi - 17
//........................................................................
nn = ijk-strideZ+strideY; // neighbor index (get convention)
m18 = Psi[nn]; // get neighbor for phi - 18
psi_Laplacian = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*phi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*psi));//Laplacian of electric potential
residual_error = psi_Laplacian+rho_e/epsilon_LB;
ResidualError[n] = residual_error;
}
}
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
// int strideY, int strideZ,int start, int finish, int Np){
//

42
models/PoissonSolver.cpp
View File

@ -338,6 +338,7 @@ void ScaLBL_Poisson::Create(){
ScaLBL_AllocateDeviceMemory((void **) &fq, 7*dist_mem_size);
ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);
ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &ResidualError, sizeof(double)*Np);
//...........................................................................
// Update GPU data structures
@ -561,31 +562,24 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
// Check convergence of steady-state solution
if (timestep%analysis_interval==0){
//ScaLBL_Comm->RegularLayout(Map,Psi,Psi_host);
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
double count_loc=0;
double count;
double psi_avg;
double psi_loc=0.f;
ScaLBL_D3Q7_PoissonResidualError(NeighborList,dvcMap,ResidualError,Psi,ChargeDensity,epsilon_LB,Nx,Nx*Ny,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior())
ScaLBL_D3Q7_PoissonResidualError(NeighborList,dvcMap,ResidualError,Psi,ChargeDensity,epsilon_LB,Nx,Nx*Ny,0, ScaLBL_Comm->LastExterior())
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){
psi_loc += Psi_host(i,j,k);
count_loc+=1.0;
}
}
}
}
psi_avg=Dm->Comm.sumReduce( psi_loc);
count=Dm->Comm.sumReduce( count_loc);
psi_avg /= count;
double psi_avg_mag=psi_avg;
if (psi_avg==0.0) psi_avg_mag=1.0;
error = fabs(psi_avg-psi_avg_previous)/fabs(psi_avg_mag);
psi_avg_previous = psi_avg;
vector<double> ResidualError_host(Np);
double error_loc_max;
//calculate the maximum residual error
ScaLBL_CopyToHost(ResidualError_host,ResidualError,sizeof(double)*Np);
vector<double>::iterator it_max = max_element(ResidualError_host[0],ResidualError_host[ScaLBL_Comm->LastExterior()]);
unsigned int idx_max1 = distance(ResidualError_host.begin(),it_max);
it_max = max_element(ResidualError_host[ScaLBL_Comm->FirstInterior()], ResidualError_host[ScaLBL_Comm->LastInterior()]);
unsigned int idx_max2 = distance(ResidualError_host.begin(),it_max);
if (ResidualError_host[idx_max1]>ResidualError_host[idx_max2]){
error_loc_max=ResidualError_host[idx_max1];
}
else{
error_loc_max=ResidualError_host[idx_max2];
}
error = Dm->Comm.maxReduce(error_loc_max);
}
}
if(WriteLog==true){

1
models/PoissonSolver.h
View File

@ -85,6 +85,7 @@ public:
double *Psi;
double *ElectricField;
double *ChargeDensityDummy;// for debugging
double *ResidualError;
private:
Utilities::MPI comm;