From 2b2bdee447e12ae30fd6358049ceb270a52820a1 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 16 Aug 2022 11:23:30 -0400 Subject: [PATCH 01/16] fix error in gpu poisson --- cuda/Poisson.cu | 1 + hip/Poisson.hip | 2 ++ models/PoissonSolver.cpp | 1 - 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index 1e3b75a1..29166692 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -559,6 +559,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, idx = Map[n]; Psi[idx] = psi; + Error[n] = error; // q = 0 dist[n] = W0*psi;// diff --git a/hip/Poisson.hip b/hip/Poisson.hip index 2808ffa1..a58d17c0 100644 --- a/hip/Poisson.hip +++ b/hip/Poisson.hip @@ -519,6 +519,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, idx = Map[n]; Psi[idx] = psi; + Error[n] = error; + // q = 0 dist[n] = W0*psi;// diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 4dab7323..664e23e1 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -808,7 +808,6 @@ 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 From b423bcb42b012e2418fec90e596106f6a0bd4115 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 21 Aug 2022 08:53:05 -0400 Subject: [PATCH 02/16] minor cleanup of gpu poisson --- cuda/Poisson.cu | 22 ++++++++++++++-------- hip/Poisson.hip | 26 ++++++++++++++++---------- 2 files changed, 30 insertions(+), 18 deletions(-) diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index 29166692..c9803271 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -3,7 +3,7 @@ //#include #define NBLOCKS 1024 -#define NTHREADS 256 +#define NTHREADS 512 __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(int *neighborList,int *Map, double *dist, double *Psi, int start, int finish, int Np){ int n; @@ -328,7 +328,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, f16, f17, f18; int nr1, nr2, nr3, nr4, nr5, nr6, nr7, nr8, nr9, nr10, nr11, nr12, nr13, nr14, nr15, nr16, nr17, nr18; - double error,sum_q; + double sum_q; double rlx = 1.0 / tau; int idx; @@ -545,6 +545,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, f17 = dist[18 * Np + n]; f18 = dist[17 * Np + n]; + /* Ex = (f1 - f2) * rlx * + 4.0; //NOTE the unit of electric field here is V/lu + Ey = (f3 - f4) * rlx * + 4.0; //factor 4.0 is D3Q7 lattice squared speed of sound + Ez = (f5 - f6) * rlx * 4.0; + */ Ex = (f1 - f2 + 0.5*(f7 - f8 + f9 - f10 + f11 - f12 + f13 - f14))*4.0; //NOTE the unit of electric field here is V/lu Ey = (f3 - f4 + 0.5*(f7 - f8 - f9 + f10 + f15 - f16 + f17 - f18))*4.0; Ez = (f5 - f6 + 0.5*(f11 - f12 - f13 + f14 + f15 - f16 - f17 + f18))*4.0; @@ -554,13 +560,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, sum_q = f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11+f12+f13+f14+f15+f16+f17+f18; error = 8.0*(sum_q - f0) + rho_e; - - psi = 2.0*(f0*(1.0 - rlx) + rlx*(sum_q + 0.125*rho_e)); - - idx = Map[n]; - Psi[idx] = psi; + Error[n] = error; + psi = 2.0*(f0*(1.0 - rlx) + rlx*(sum_q + 0.125*rho_e)); + + idx = Map[n]; + Psi[idx] = psi; + // q = 0 dist[n] = W0*psi;// @@ -594,7 +601,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, dist[16 * Np + n] = W2*psi;//f16 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; dist[17 * Np + n] = W2*psi;//f17 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; dist[18 * Np + n] = W2*psi;//f18 * (1.0 - rlx) +W2* (rlx * psi) - (1.0-0.5*rlx)*0.02777777777777778*rho_e; - //........................................................................ } } diff --git a/hip/Poisson.hip b/hip/Poisson.hip index a58d17c0..01844b60 100644 --- a/hip/Poisson.hip +++ b/hip/Poisson.hip @@ -301,7 +301,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, //........Get 1-D index for this thread.................... n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start; if (n Date: Mon, 22 Aug 2022 22:00:19 -0400 Subject: [PATCH 03/16] fix cuda thing --- cuda/Poisson.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index c9803271..fb8f716a 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -421,7 +421,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, f18 = dist[nr18]; sum_q = f1+f2+f3+f4+f5+f6+f7+f8+f9+f10+f11+f12+f13+f14+f15+f16+f17+f18; - error = 8.0*(sum_q - f0) + rho_e; + //error = 8.0*(sum_q - f0) + rho_e; psi = 2.0*(f0*(1.0 - rlx) + rlx*(sum_q + 0.125*rho_e)); From 372ebd3af201dcb25cb7fa0804cc09c2bb95d3c4 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 25 Aug 2022 15:46:47 -0400 Subject: [PATCH 04/16] add flat field to accelerate Poisson --- models/PoissonSolver.cpp | 139 ++++++++++++++++++++ models/PoissonSolver.h | 2 + tests/lbpm_nernst_planck_cell_simulator.cpp | 2 +- 3 files changed, 142 insertions(+), 1 deletion(-) diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 664e23e1..4f611204 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -862,6 +862,145 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times } } +void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bool UseSlippingVelBC, int timestep_from_Study){ + + double error = 1.0; + double threshold = 10000000.0; + bool SET_THRESHOLD = false; + if (electric_db->keyExists( "rescale_at_distance" )){ + SET_THRESHOLD = true; + threshold = electric_db->getScalar( "rescale_at_distance" ); + } + if (BoundaryConditionInlet > 0) SET_THRESHOLD = false; + if (BoundaryConditionOutlet > 0) SET_THRESHOLD = false; + + double *host_Error; + host_Error = new double [Np]; + + timestep=0; + auto t1 = std::chrono::system_clock::now(); + while (timestep < timestepMax && error > tolerance) { + //************************************************************************/ + // *************ODD TIMESTEP*************// + timestep++; + //SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential + SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision + ScaLBL_Comm->Barrier(); comm.barrier(); + + // *************EVEN TIMESTEP*************// + timestep++; + //SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential + SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision + ScaLBL_Comm->Barrier(); comm.barrier(); + //************************************************************************/ + + // Check convergence of steady-state solution + if (timestep==2){ + //save electric potential for convergence check + } + if (timestep%analysis_interval==0){ + /* get the elecric potential */ + ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz); + if (rank==0) printf(" ... getting Poisson solver error \n"); + double err = 0.0; + double max_error = 0.0; + ScaLBL_CopyToHost(host_Error,ResidualError,sizeof(double)*Np); + for (int idx=0; idx max_error ){ + max_error = err; + } + } + error=Dm->Comm.maxReduce(max_error); + + if (error > tolerance & SET_THRESHOLD){ + /* don't use this with an external BC */ + // cpompute the far-field electric potential + double inside_local = 0.0; + double outside_local = 0.0; + double inside_count_local = 0.0; + double outside_count_local = 0.0; + /* global values */ + double inside_global = 0.0; + double outside_global = 0.0; + double inside_count_global = 0.0; + double outside_count_global = 0.0; + for (int k=1; k threshold && distance < (threshold + 1.0)){ + outside_count_local += 1.0; + outside_local += Psi_host(n); + } + else if (distance < (-1.0)*threshold && distance > (-1.0)*(threshold + 1.0)){ + inside_count_local += 1.0; + inside_local += Psi_host(n); + } + } + } + } + inside_count_global = Dm->Comm.sumReduce(inside_count_local); + outside_count_global = Dm->Comm.sumReduce(outside_count_local); + outside_global = Dm->Comm.sumReduce(outside_local); + inside_global = Dm->Comm.sumReduce(inside_local); + outside_global /= outside_count_global; + inside_global /= inside_count_global; + + if (rank==0) printf(" Rescaling far-field electric potential to value (outside): %f \n",outside_global); + if (rank==0) printf(" Rescaling far-field electric potential to value (inside): %f \n",inside_global); + + // rescale the far-field electric potential + for (int k=1; k (threshold + 1.0)){ + Psi_host(n) = outside_global; + } + else if ( distance < (-1.0)*(threshold + 1.0)){ + Psi_host(n) = inside_global; + } + } + } + } + ScaLBL_CopyToDevice(Psi,Psi_host.data(),sizeof(double)*Nx*Ny*Nz); + } + /* compute the eletric field */ + //ScaLBL_D3Q19_Poisson_getElectricField(fq, ElectricField, tau, Np); + } + } + if (rank == 0) + printf("---------------------------------------------------------------" + "----\n"); + // Compute the walltime per timestep + auto t2 = std::chrono::system_clock::now(); + double cputime = std::chrono::duration(t2 - t1).count() / timestep; + // Performance obtained from each node + double MLUPS = double(Np) / cputime / 1000000; + + if (rank == 0) + printf("********************************************************\n"); + if (rank == 0) + printf("CPU time = %f \n", cputime); + if (rank == 0) + printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; + if (rank == 0) + printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + if (rank == 0) + printf("********************************************************\n"); + delete [] host_Error; + + //************************************************************************/ + + if(WriteLog==true){ + getConvergenceLog(timestep,error); + } +} + void ScaLBL_Poisson::getConvergenceLog(int timestep,double error){ if ( rank == 0 ) { fprintf(TIMELOG,"%i %.5g\n",timestep,error); diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 8ab60140..faee8c01 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -35,6 +35,8 @@ public: void Create(); void Initialize(double time_conv_from_Study); void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study); + void Run(double *ChargeDensity, DoubleArray MembraneDistance, + bool UseSlippingVelBC, int timestep_from_Study); void getElectricPotential(DoubleArray &ReturnValues); void getElectricPotential_debug(int timestep); void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y, diff --git a/tests/lbpm_nernst_planck_cell_simulator.cpp b/tests/lbpm_nernst_planck_cell_simulator.cpp index 3ac2065c..f98ed1b7 100644 --- a/tests/lbpm_nernst_planck_cell_simulator.cpp +++ b/tests/lbpm_nernst_planck_cell_simulator.cpp @@ -125,7 +125,7 @@ int main(int argc, char **argv) while (timestep < Study.timestepMax){ timestep++; - PoissonSolver.Run(IonModel.ChargeDensity,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental + PoissonSolver.Run(IonModel.ChargeDensity,IonModel.MembraneDistance,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental //comm.barrier(); //if (rank == 0) printf(" Poisson step %i \n",timestep); //StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity From fc4af2d712e58c4594e0fcd0727d099778dbf1ed Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 27 Aug 2022 15:40:32 -0400 Subject: [PATCH 05/16] fix swc thing --- common/Domain.cpp | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 945a6c6a..14dc17dd 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -217,17 +217,23 @@ void Domain::read_swc(const std::string &Filename) { double z = k*voxel_length + start_z; double distance; - double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma); + double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma); + + double di = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi)); + double dp = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp)); + if (s > length){ - distance = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi)); + distance = di; } else if (s < 0.0){ - distance = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp)); + distance = dp; } else { // linear variation for radius double radius = rp + (ri - rp)*s/length; distance = radius - sqrt((x-xp-alpha*s)*(x-xp-alpha*s) + (y-yp-beta*s)*(y-yp-beta*s) + (z-zp-gamma*s)*(z-zp-gamma*s)); + if (distance < di) distance = di; + if (distance < dp) distance = dp; } if ( distance > 0.0 ){ /* label the voxel */ From 9229de287510fdb42af7b29939152415f308d8d9 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 27 Aug 2022 16:06:30 -0400 Subject: [PATCH 06/16] fix swc stuff try 2 --- common/Domain.cpp | 22 +++++++-------------- tests/lbpm_nernst_planck_cell_simulator.cpp | 2 +- 2 files changed, 8 insertions(+), 16 deletions(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 14dc17dd..94d41a7b 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -215,26 +215,18 @@ void Domain::read_swc(const std::string &Filename) { double x = i*voxel_length + start_x; double y = j*voxel_length + start_y; double z = k*voxel_length + start_z; - + double distance; double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma); - + double di = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi)); double dp = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp)); + // linear variation for radius + double radius = rp + (ri - rp)*s/length; + distance = radius - sqrt((x-xp-alpha*s)*(x-xp-alpha*s) + (y-yp-beta*s)*(y-yp-beta*s) + (z-zp-gamma*s)*(z-zp-gamma*s)); + if (distance < di) distance = di; + if (distance < dp) distance = dp; - if (s > length){ - distance = di; - } - else if (s < 0.0){ - distance = dp; - } - else { - // linear variation for radius - double radius = rp + (ri - rp)*s/length; - distance = radius - sqrt((x-xp-alpha*s)*(x-xp-alpha*s) + (y-yp-beta*s)*(y-yp-beta*s) + (z-zp-gamma*s)*(z-zp-gamma*s)); - if (distance < di) distance = di; - if (distance < dp) distance = dp; - } if ( distance > 0.0 ){ /* label the voxel */ //id[k*Nx*Ny + j*Nx + i] = label; diff --git a/tests/lbpm_nernst_planck_cell_simulator.cpp b/tests/lbpm_nernst_planck_cell_simulator.cpp index f98ed1b7..d038b782 100644 --- a/tests/lbpm_nernst_planck_cell_simulator.cpp +++ b/tests/lbpm_nernst_planck_cell_simulator.cpp @@ -133,7 +133,7 @@ int main(int argc, char **argv) IonModel.RunMembrane(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane //comm.barrier(); - if (rank == 0) printf(" Membrane step %i \n",timestep); + //if (rank == 0) printf(" Membrane step %i \n",timestep); fflush(stdout); From ffe4f794debbdccd2723a37effc695e1a9736845 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 27 Aug 2022 16:30:03 -0400 Subject: [PATCH 07/16] fix swc stuff try 3 --- common/Domain.cpp | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 94d41a7b..d5411464 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -215,18 +215,27 @@ void Domain::read_swc(const std::string &Filename) { double x = i*voxel_length + start_x; double y = j*voxel_length + start_y; double z = k*voxel_length + start_z; - + double distance; double s = ((x-xp)*alpha+(y-yp)*beta+(z-zp)*gamma) / (alpha*alpha + beta*beta + gamma*gamma); - + double di = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi)); double dp = rp - sqrt((x-xp)*(x-xp) + (y-yp)*(y-yp) + (z-zp)*(z-zp)); - // linear variation for radius - double radius = rp + (ri - rp)*s/length; - distance = radius - sqrt((x-xp-alpha*s)*(x-xp-alpha*s) + (y-yp-beta*s)*(y-yp-beta*s) + (z-zp-gamma*s)*(z-zp-gamma*s)); + + if (s > length ){ + distance = di; + } + else if (s < 0.0){ + distance = dp; + } + else { + // linear variation for radius + double radius = rp + (ri - rp)*s/length; + distance = radius - sqrt((x-xp-alpha*s)*(x-xp-alpha*s) + (y-yp-beta*s)*(y-yp-beta*s) + (z-zp-gamma*s)*(z-zp-gamma*s)); + } if (distance < di) distance = di; if (distance < dp) distance = dp; - + if ( distance > 0.0 ){ /* label the voxel */ //id[k*Nx*Ny + j*Nx + i] = label; From e9096dbfc3c003f8e8075b46b6d8525c031a61b8 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 31 Aug 2022 14:08:59 -0400 Subject: [PATCH 08/16] update docs --- docs/source/userGuide/models/color/index.rst | 40 +++++++++++++++--- .../models/greyscaleColor/greyscaleColor.rst | 41 ++++++++++++++++--- 2 files changed, 69 insertions(+), 12 deletions(-) diff --git a/docs/source/userGuide/models/color/index.rst b/docs/source/userGuide/models/color/index.rst index 8a7f2519..47ca36a7 100644 --- a/docs/source/userGuide/models/color/index.rst +++ b/docs/source/userGuide/models/color/index.rst @@ -176,42 +176,70 @@ The non-zero equilibrium moments are defined as :nowrap: $$ - m_1^{eq} = (j_x^2+j_y^2+j_z^2) - \alpha |\textbf{C}|, \\ + m_1^{eq} = 19\frac{ j_x^2+j_y^2+j_z^2}{\rho_0} - 11\rho - 19 \alpha |\textbf{C}|, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_2^{eq} = 3\rho - \frac{11( j_x^2+j_y^2+j_z^2)}{2\rho_0}, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_4^{eq} = -\frac{2 j_x}{3}, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_6^{eq} = -\frac{2 j_y}{3}, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_8^{eq} = -\frac{2 j_z}{3}, \\ $$ .. math:: :nowrap: $$ - m_9^{eq} = (2j_x^2-j_y^2-j_z^2)+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\ + m_9^{eq} = \frac{2j_x^2-j_y^2-j_z^2}{\rho_0}+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\ $$ .. math:: :nowrap: $$ - m_{11}^{eq} = (j_y^2-j_z^2) + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\ + m_{11}^{eq} = \frac{j_y^2-j_z^2}{\rho_0} + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\ $$ .. math:: :nowrap: $$ - m_{13}^{eq} = j_x j_y + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\ + m_{13}^{eq} = \frac{j_x j_y}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\ $$ .. math:: :nowrap: $$ - m_{14}^{eq} = j_y j_z + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\ + m_{14}^{eq} = \frac{j_y j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\ $$ .. math:: :nowrap: $$ - m_{15}^{eq} = j_x j_z + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;, + m_{15}^{eq} = \frac{j_x j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;. $$ where the color gradient is determined from the phase indicator field diff --git a/docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst b/docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst index a86aa5fc..46094d75 100644 --- a/docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst +++ b/docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst @@ -241,46 +241,75 @@ The relaxation parameters are determined from the relaxation time: The non-zero equilibrium moments are defined as + .. math:: :nowrap: $$ - m_1^{eq} = (j_x^2+j_y^2+j_z^2) - \alpha |\textbf{C}|, \\ + m_1^{eq} = 19\frac{ j_x^2+j_y^2+j_z^2}{\rho_0} - 11\rho - 19 \alpha |\textbf{C}|, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_2^{eq} = 3\rho - \frac{11( j_x^2+j_y^2+j_z^2)}{2\rho_0}, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_4^{eq} = -\frac{2 j_x}{3}, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_6^{eq} = -\frac{2 j_y}{3}, \\ + $$ + +.. math:: + :nowrap: + + $$ + m_8^{eq} = -\frac{2 j_z}{3}, \\ $$ .. math:: :nowrap: $$ - m_9^{eq} = (2j_x^2-j_y^2-j_z^2)+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\ + m_9^{eq} = \frac{2j_x^2-j_y^2-j_z^2}{\rho_0}+ \alpha \frac{|\textbf{C}|}{2}(2n_x^2-n_y^2-n_z^2), \\ $$ .. math:: :nowrap: $$ - m_{11}^{eq} = (j_y^2-j_z^2) + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\ + m_{11}^{eq} = \frac{j_y^2-j_z^2}{\rho_0} + \alpha \frac{|\textbf{C}|}{2}(n_y^2-n_z^2), \\ $$ .. math:: :nowrap: $$ - m_{13}^{eq} = j_x j_y + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\ + m_{13}^{eq} = \frac{j_x j_y}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_y\;, \\ $$ .. math:: :nowrap: $$ - m_{14}^{eq} = j_y j_z + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\ + m_{14}^{eq} = \frac{j_y j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_y n_z\;, \\ $$ .. math:: :nowrap: $$ - m_{15}^{eq} = j_x j_z + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;, + m_{15}^{eq} = \frac{j_x j_z}{\rho_0} + \alpha \frac{|\textbf{C}|}{2} n_x n_z\;. $$ where the color gradient is determined from the phase indicator field From f1ceb930ee4acc8710d2fb8206b0ed34a4e9260c Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 26 Oct 2022 16:16:59 -0400 Subject: [PATCH 09/16] working on read/write for membrane coefficients --- common/Membrane.cpp | 68 +++++++++++++++++++++++++++++++++++++++++++++ common/Membrane.h | 12 ++++++++ 2 files changed, 80 insertions(+) diff --git a/common/Membrane.cpp b/common/Membrane.cpp index ccc530cf..c9943b42 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -667,6 +667,74 @@ int Membrane::Create(DoubleArray &Distance, IntArray &Map){ return mlink; } +void Membrane::Write(string filename){ + + int mlink = membraneLinkCount; + std::ofstream ofs (filename, std::ofstream::out); + /* Create local copies of membrane data structures */ + double *tmpMembraneCoef; // mass transport coefficient for the membrane + tmpMembraneCoef = new double [2*mlink*sizeof(double)]; + ScaLBL_CopyToHost(tmpMembraneCoef, MembraneCoef, 2*mlink*sizeof(double)); + int i,j,k; + for (int m=0; m Date: Wed, 26 Oct 2022 16:17:23 -0400 Subject: [PATCH 10/16] 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( "tolerance_method", "MSE" ); - lattice_scheme = electric_db->getWithDefault( "lattice_scheme", "D3Q7" ); + lattice_scheme = electric_db->getWithDefault( "lattice_scheme", "D3Q19" ); if (electric_db->keyExists( "epsilonR" )){ epsilonR = electric_db->getScalar( "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( "format", "hdf5" ); + + DoubleArray ElectricalPotential(Nx, Ny, Nz); + std::vector visData; + fillHalo 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(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(); + + //-------------------------------------------------------------------------------------------------------------------- + + //-------------------------------------Create Names for Variables------------------------------------------------------ + if (vis_db->getWithDefault("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("save_electric_potential", true)) { + ASSERT(visData[0].vars[0]->name == "ElectricPotential"); + getElectricPotential(ElectricalPotential); + Array &ElectricPotentialData = visData[0].vars[0]->data; + fillData.copy(ElectricalPotential, ElectricPotentialData); + } + + if (vis_db->getWithDefault("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 From 79484ab8fa041ce0cce8e3fd2f8c67a9bdea0f1e Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 26 Oct 2022 20:24:45 -0400 Subject: [PATCH 11/16] update poisson solver --- tests/TestPoissonSolver.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/TestPoissonSolver.cpp b/tests/TestPoissonSolver.cpp index 7d39889a..637755d0 100644 --- a/tests/TestPoissonSolver.cpp +++ b/tests/TestPoissonSolver.cpp @@ -76,11 +76,14 @@ int main(int argc, char **argv) PoissonSolver.getElectricField_debug(timestep); } } + PoissonSolver.WriteVis(timestep); } else { + int timestep = 1; PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,1); PoissonSolver.getElectricPotential_debug(1); PoissonSolver.getElectricField_debug(1); + PoissonSolver.WriteVis(timestep); } if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); From 95233d61a7a00ad15b324474ebba584ff5d6eef8 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 26 Oct 2022 20:24:56 -0400 Subject: [PATCH 12/16] update configure script --- sample_scripts/configure_arden | 1 + 1 file changed, 1 insertion(+) diff --git a/sample_scripts/configure_arden b/sample_scripts/configure_arden index 58d303f9..988300cd 100755 --- a/sample_scripts/configure_arden +++ b/sample_scripts/configure_arden @@ -13,6 +13,7 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc \ -D CUDA_HOST_COMPILER="/usr/bin/gcc" \ -D USE_HDF5=1 \ -D HDF5_DIRECTORY="/opt/arden/hdf5/1.8.12" \ + -D USE_DOXYGEN=true \ -D USE_CUDA=0 \ -D USE_TIMER=0 \ ~/Programs/LBPM-WIA From 111529ff5ef485497b4f719fda7d459dd705dc27 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 27 Oct 2022 08:15:21 -0400 Subject: [PATCH 13/16] add voltage BC for d3q19 poisson solver --- common/ScaLBL.cpp | 23 +++++++ common/ScaLBL.h | 8 +++ cpu/D3Q7BC.cpp | 3 + cpu/Poisson.cpp | 88 ++++++++++++++++++++++++++ models/IonModel.cpp | 3 +- models/PoissonSolver.cpp | 72 +++++++++++++++++---- models/PoissonSolver.h | 4 +- tests/lbpm_nernst_planck_simulator.cpp | 1 - 8 files changed, 186 insertions(+), 16 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index f83bcb85..a865b4aa 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -2472,6 +2472,29 @@ void ScaLBL_Communicator::D3Q7_Poisson_Potential_BC_Z(int *neighborList, double } } + +void ScaLBL_Communicator::D3Q19_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time){ + if (kproc == 0) { + if (time%2==0){ + ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(dvcSendList_z, fq, Vin, sendCount_z, N); + } + else{ + ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(neighborList, dvcSendList_z, fq, Vin, sendCount_z, N); + } + } +} + +void ScaLBL_Communicator::D3Q19_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time){ + if (kproc == nprocz-1){ + if (time%2==0){ + ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(dvcSendList_Z, fq, Vout, sendCount_Z, N); + } + else{ + ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(neighborList, dvcSendList_Z, fq, Vout, sendCount_Z, N); + } + } +} + void ScaLBL_Communicator::Poisson_D3Q7_BC_z(int *Map, double *Psi, double Vin){ if (kproc == 0) { ScaLBL_Poisson_D3Q7_BC_z(dvcSendList_z, Map, Psi, Vin, sendCount_z); diff --git a/common/ScaLBL.h b/common/ScaLBL.h index cc69302e..64d178d5 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -396,6 +396,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, extern "C" void ScaLBL_D3Q19_Poisson_getElectricField(double *dist, double *ElectricField, double tau, int Np); + +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np); +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np); +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np); +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vout, int count, int Np); + // LBM Stokes Model (adapted from MRT model) extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np); @@ -806,6 +812,8 @@ public: double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time); void D3Q7_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time); void D3Q7_Poisson_Potential_BC_Z(int *neighborList, double *fq, double Vout, int time); + void D3Q19_Poisson_Potential_BC_z(int *neighborList, double *fq, double Vin, int time); + void D3Q19_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); diff --git a/cpu/D3Q7BC.cpp b/cpu/D3Q7BC.cpp index c99bd257..3f30fde4 100644 --- a/cpu/D3Q7BC.cpp +++ b/cpu/D3Q7BC.cpp @@ -181,6 +181,7 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_z(int *d_neighborList, } } + extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, @@ -213,6 +214,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, } } + + extern "C" void ScaLBL_Poisson_D3Q7_BC_z(int *list, int *Map, double *Psi, double Vin, int count) { int idx, n, nm; diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp index 9f71bb6d..97ee20e7 100644 --- a/cpu/Poisson.cpp +++ b/cpu/Poisson.cpp @@ -631,6 +631,8 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_ElectricPotential( } } + + extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, @@ -915,6 +917,92 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson(int *Map, double *dist, } } +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) { + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int nread, nr5; + + for (int idx = 0; idx < count; idx++) { + int n = list[idx]; + + dist[6 * Np + n] = W1*Vin; + dist[12 * Np + n] = W2*Vin; + dist[13 * Np + n] = W2*Vin; + dist[16 * Np + n] = W2*Vin; + dist[17 * Np + n] = W2*Vin; + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, + double *dist, + double Vout, + int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + + for (int idx = 0; idx < count; idx++) { + + int n = list[idx]; + dist[5 * Np + n] = W1*Vout; + dist[11 * Np + n] = W2*Vout; + dist[14 * Np + n] = W2*Vout; + dist[15 * Np + n] = W2*Vout; + dist[18 * Np + n] = W2*Vout; + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, + int *list, + double *dist, + double Vin, int count, + int Np) { + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int nr5, nr11, nr14, nr15, nr18; + + for (int idx = 0; idx < count; idx++) { + int n = list[idx]; + + // Unknown distributions + nr5 = d_neighborList[n + 4 * Np]; + nr11 = d_neighborList[n + 10 * Np]; + nr15 = d_neighborList[n + 14 * Np]; + nr14 = d_neighborList[n + 13 * Np]; + nr18 = d_neighborList[n + 17 * Np]; + + dist[nr5] = W1*Vin; + dist[nr11] = W2*Vin; + dist[nr15] = W2*Vin; + dist[nr14] = W2*Vin; + dist[nr18] = W2*Vin; + } +} + +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int nr6, nr12, nr13, nr16, nr17; + + for (int idx = 0; idx < count; idx++) { + + int n = list[idx]; + // unknown distributions + nr6 = d_neighborList[n + 5 * Np]; + nr12 = d_neighborList[n + 11 * Np]; + nr16 = d_neighborList[n + 15 * Np]; + nr17 = d_neighborList[n + 16 * Np]; + nr13 = d_neighborList[n + 12 * Np]; + + dist[nr6] = W1*Vout; + dist[nr12] = W2*Vout; + dist[nr16] = W2*Vout; + dist[nr17] = W2*Vout; + dist[nr13] = W2*Vout; + } +} + extern "C" void ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np) { int n; diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 4e389de3..16679884 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -300,7 +300,8 @@ void ScaLBL_IonModel::ReadParams(string filename, vector &num_iter) { 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 - USE_MEMBRANE = true; + USE_MEMBRANE = false; + Restart = false; // read the input database db = std::make_shared(filename); domain_db = db->getDatabase("Domain"); diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 89bfb985..ee4ccf1e 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -726,13 +726,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times // *************ODD TIMESTEP*************// timestep++; SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential - SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision + SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); // *************EVEN TIMESTEP*************// timestep++; SolveElectricPotentialAAeven(timestep_from_Study);//update electric potential - SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision + SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); //************************************************************************/ @@ -797,14 +797,14 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times //************************************************************************/ // *************ODD TIMESTEP*************// timestep++; - //SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision + //SolveElectricPotentialAAodd(timestep_from_Study);//,ChargeDensity, UseSlippingVelBC);//update electric potential + SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); // *************EVEN TIMESTEP*************// timestep++; - //SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision + //SolveElectricPotentialAAeven(timestep_from_Study);//,ChargeDensity, UseSlippingVelBC);//update electric potential + SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); //************************************************************************/ @@ -884,13 +884,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bo // *************ODD TIMESTEP*************// timestep++; //SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision + SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC,timestep);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); // *************EVEN TIMESTEP*************// timestep++; //SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential - SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision + SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); //************************************************************************/ @@ -913,7 +913,7 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bo } error=Dm->Comm.maxReduce(max_error); - if (error > tolerance & SET_THRESHOLD){ + if (error > tolerance && SET_THRESHOLD){ /* don't use this with an external BC */ // cpompute the far-field electric potential double inside_local = 0.0; @@ -1146,7 +1146,7 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ } } -void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC){ +void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC, int timestep){ if (lattice_scheme.compare("D3Q7")==0){ ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -1167,6 +1167,30 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVe ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); //ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + + // Set boundary conditions + if (BoundaryConditionInlet > 0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } + } ScaLBL_D3Q19_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); @@ -1176,7 +1200,7 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVe } } -void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC){ +void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC, int timestep){ if (lattice_scheme.compare("D3Q7")==0){ ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -1194,6 +1218,31 @@ void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingV ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE // ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + + // Set boundary conditions + if (BoundaryConditionInlet > 0){ + switch (BoundaryConditionInlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + case 2: + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); + break; + } + } + if (BoundaryConditionOutlet > 0){ + switch (BoundaryConditionOutlet){ + case 1: + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + case 2: + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep); + ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); + break; + } + } + ScaLBL_D3Q19_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, ResidualError, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); @@ -1323,7 +1372,6 @@ 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( "format", "hdf5" ); DoubleArray ElectricalPotential(Nx, Ny, Nz); diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index c6672cec..e712f76e 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -118,8 +118,8 @@ private: void SolveElectricPotentialAAeven(int timestep_from_Study); void SolveElectricPotentialAAodd(int timestep_from_Study); //void SolveElectricField(); - void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC); - void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC); + void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC, int timestep); + void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC, int timestep); void getConvergenceLog(int timestep,double error); double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step); diff --git a/tests/lbpm_nernst_planck_simulator.cpp b/tests/lbpm_nernst_planck_simulator.cpp index bb2389cd..033adfd3 100644 --- a/tests/lbpm_nernst_planck_simulator.cpp +++ b/tests/lbpm_nernst_planck_simulator.cpp @@ -104,7 +104,6 @@ int main(int argc, char **argv) //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); //} From 14c1037d7da43ef18c80db567f9510ee3e708773 Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 28 Oct 2022 07:20:45 -0400 Subject: [PATCH 14/16] add GPU voltage bc --- cuda/Poisson.cu | 125 ++++++++++++++++++++++++++++++++++++++++++++++++ hip/Poisson.hip | 125 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 250 insertions(+) diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index fb8f716a..04719551 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -642,6 +642,131 @@ __global__ void dvc_ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *P } } +__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + + dist[6 * Np + n] = W1*Vin; + dist[12 * Np + n] = W2*Vin; + dist[13 * Np + n] = W2*Vin; + dist[16 * Np + n] = W2*Vin; + dist[17 * Np + n] = W2*Vin; + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + dist[5 * Np + n] = W1*Vout; + dist[11 * Np + n] = W2*Vout; + dist[14 * Np + n] = W2*Vout; + dist[15 * Np + n] = W2*Vout; + dist[18 * Np + n] = W2*Vout; + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int nr5, nr11, nr14, nr15, nr18; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + } + + // Unknown distributions + nr5 = d_neighborList[n + 4 * Np]; + nr11 = d_neighborList[n + 10 * Np]; + nr15 = d_neighborList[n + 14 * Np]; + nr14 = d_neighborList[n + 13 * Np]; + nr18 = d_neighborList[n + 17 * Np]; + + dist[nr5] = W1*Vin; + dist[nr11] = W2*Vin; + dist[nr15] = W2*Vin; + dist[nr14] = W2*Vin; + dist[nr18] = W2*Vin; + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int nr6, nr12, nr13, nr16, nr17; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + // unknown distributions + nr6 = d_neighborList[n + 5 * Np]; + nr12 = d_neighborList[n + 11 * Np]; + nr16 = d_neighborList[n + 15 * Np]; + nr17 = d_neighborList[n + 16 * Np]; + nr13 = d_neighborList[n + 12 * Np]; + + dist[nr6] = W1*Vout; + dist[nr12] = W2*Vout; + dist[nr16] = W2*Vout; + dist[nr17] = W2*Vout; + dist[nr13] = W2*Vout; + } +} + +/* wrapper functions to launch kernels */ +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vout, count, N); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err)); + } +} +// +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<>>(list, dist, Vout, count, N); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err)); + } +} +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count,int Np) { + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<>>(d_neighborList, list, dist, Vout, count, N); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err)); + } +} +// +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) { + + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<>>(d_neighborList, list, dist, Vout, count, N); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err)); + } +} + + extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, diff --git a/hip/Poisson.hip b/hip/Poisson.hip index 01844b60..d1c666a4 100644 --- a/hip/Poisson.hip +++ b/hip/Poisson.hip @@ -603,6 +603,131 @@ __global__ void dvc_ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *P } } +__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + + dist[6 * Np + n] = W1*Vin; + dist[12 * Np + n] = W2*Vin; + dist[13 * Np + n] = W2*Vin; + dist[16 * Np + n] = W2*Vin; + dist[17 * Np + n] = W2*Vin; + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + dist[5 * Np + n] = W1*Vout; + dist[11 * Np + n] = W2*Vout; + dist[14 * Np + n] = W2*Vout; + dist[15 * Np + n] = W2*Vout; + dist[18 * Np + n] = W2*Vout; + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int nr5, nr11, nr14, nr15, nr18; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + } + + // Unknown distributions + nr5 = d_neighborList[n + 4 * Np]; + nr11 = d_neighborList[n + 10 * Np]; + nr15 = d_neighborList[n + 14 * Np]; + nr14 = d_neighborList[n + 13 * Np]; + nr18 = d_neighborList[n + 17 * Np]; + + dist[nr5] = W1*Vin; + dist[nr11] = W2*Vin; + dist[nr15] = W2*Vin; + dist[nr14] = W2*Vin; + dist[nr18] = W2*Vin; + } +} + +__global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) { + + double W1 = 1.0/24.0; + double W2 = 1.0/48.0; + int nr6, nr12, nr13, nr16, nr17; + + int idx = blockIdx.x*blockDim.x + threadIdx.x; + + if (idx < count){ + int n = list[idx]; + // unknown distributions + nr6 = d_neighborList[n + 5 * Np]; + nr12 = d_neighborList[n + 11 * Np]; + nr16 = d_neighborList[n + 15 * Np]; + nr17 = d_neighborList[n + 16 * Np]; + nr13 = d_neighborList[n + 12 * Np]; + + dist[nr6] = W1*Vout; + dist[nr12] = W2*Vout; + dist[nr16] = W2*Vout; + dist[nr17] = W2*Vout; + dist[nr13] = W2*Vout; + } +} + +/* wrapper functions to launch kernels */ +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vout, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("HIP error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z (kernel): %s \n",hipGetErrorString(err)); + } +} +// +extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){ + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<>>(list, dist, Vout, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("HIP error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z (kernel): %s \n",hipGetErrorString(err)); + } +} +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count,int Np) { + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<>>(d_neighborList, list, dist, Vout, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("HIP error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z (kernel): %s \n",hipGetErrorString(err)); + } +} +// +extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) { + + int GRID = count / 512 + 1; + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<>>(d_neighborList, list, dist, Vout, count, N); + hipError_t err = hipGetLastError(); + if (hipSuccess != err){ + printf("HIP error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z (kernel): %s \n",hipGetErrorString(err)); + } +} + + extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, From 9749600540c482c569c87e4646de9c460edaa1dd Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 29 Oct 2022 13:55:30 -0400 Subject: [PATCH 15/16] fix poisson cuda --- cuda/Poisson.cu | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index 04719551..c86043a5 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -681,13 +681,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborLi double W1 = 1.0/24.0; double W2 = 1.0/48.0; - int nr5, nr11, nr14, nr15, nr18; + int nr5, nr11, nr14, nr15, nr18; int idx = blockIdx.x*blockDim.x + threadIdx.x; if (idx < count){ int n = list[idx]; - } + // Unknown distributions nr5 = d_neighborList[n + 4 * Np]; @@ -708,7 +708,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborLi double W1 = 1.0/24.0; double W2 = 1.0/48.0; - int nr6, nr12, nr13, nr16, nr17; + int nr6, nr12, nr13, nr16, nr17; int idx = blockIdx.x*blockDim.x + threadIdx.x; @@ -732,7 +732,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborLi /* wrapper functions to launch kernels */ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){ int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vout, count, N); + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vin, count, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err)); @@ -741,7 +741,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *d // extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){ int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<>>(list, dist, Vout, count, N); + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<>>(list, dist, Vout, count, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err)); @@ -749,7 +749,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *di } extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count,int Np) { int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<>>(d_neighborList, list, dist, Vout, count, N); + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<>>(d_neighborList, list, dist, Vin, count, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err)); @@ -759,7 +759,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, i extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) { int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<>>(d_neighborList, list, dist, Vout, count, N); + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<>>(d_neighborList, list, dist, Vout, count, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err)); From 17698b7dd5a240f4dcad81081c8d9f0d4f086d51 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 29 Oct 2022 13:57:54 -0400 Subject: [PATCH 16/16] fix for hip poisson --- hip/Poisson.hip | 40 +++++++++++++++++++--------------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/hip/Poisson.hip b/hip/Poisson.hip index d1c666a4..ac36000d 100644 --- a/hip/Poisson.hip +++ b/hip/Poisson.hip @@ -602,7 +602,6 @@ __global__ void dvc_ScaLBL_D3Q19_Poisson_Init(int *Map, double *dist, double *P } } } - __global__ void dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) { double W1 = 1.0/24.0; @@ -642,13 +641,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborLi double W1 = 1.0/24.0; double W2 = 1.0/48.0; - int nr5, nr11, nr14, nr15, nr18; + int nr5, nr11, nr14, nr15, nr18; int idx = blockIdx.x*blockDim.x + threadIdx.x; if (idx < count){ int n = list[idx]; - } + // Unknown distributions nr5 = d_neighborList[n + 4 * Np]; @@ -669,7 +668,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborLi double W1 = 1.0/24.0; double W2 = 1.0/48.0; - int nr6, nr12, nr13, nr16, nr17; + int nr6, nr12, nr13, nr16, nr17; int idx = blockIdx.x*blockDim.x + threadIdx.x; @@ -693,41 +692,40 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborLi /* wrapper functions to launch kernels */ extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){ int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vout, count, N); - hipError_t err = hipGetLastError(); - if (hipSuccess != err){ - printf("HIP error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z (kernel): %s \n",hipGetErrorString(err)); + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vin, count, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("Hip error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err)); } } // extern "C" void ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z(int *list, double *dist, double Vout, int count, int Np){ int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<>>(list, dist, Vout, count, N); - hipError_t err = hipGetLastError(); - if (hipSuccess != err){ - printf("HIP error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z (kernel): %s \n",hipGetErrorString(err)); + dvc_ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z<<>>(list, dist, Vout, count, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("Hip error in ScaLBL_D3Q19_AAeven_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err)); } } extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z(int *d_neighborList, int *list, double *dist, double Vin, int count,int Np) { int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<>>(d_neighborList, list, dist, Vout, count, N); - hipError_t err = hipGetLastError(); - if (hipSuccess != err){ - printf("HIP error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z (kernel): %s \n",hipGetErrorString(err)); + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z<<>>(d_neighborList, list, dist, Vin, count, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("Hip error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_z (kernel): %s \n",cudaGetErrorString(err)); } } // extern "C" void ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z(int *d_neighborList, int *list, double *dist, double Vout, int count, int Np) { int GRID = count / 512 + 1; - dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<>>(d_neighborList, list, dist, Vout, count, N); - hipError_t err = hipGetLastError(); - if (hipSuccess != err){ - printf("HIP error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z (kernel): %s \n",hipGetErrorString(err)); + dvc_ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z<<>>(d_neighborList, list, dist, Vout, count, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("Hip error in ScaLBL_D3Q19_AAodd_Poisson_Potential_BC_Z (kernel): %s \n",cudaGetErrorString(err)); } } - extern "C" void ScaLBL_D3Q19_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField,