diff --git a/common/Domain.cpp b/common/Domain.cpp index 945a6c6a..d5411464 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -217,18 +217,25 @@ 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); - if (s > length){ - distance = ri - sqrt((x-xi)*(x-xi) + (y-yi)*(y-yi) + (z-zi)*(z-zi)); + 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 = 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 */ //id[k*Nx*Ny + j*Nx + i] = label; 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 #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; @@ -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)); @@ -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,12 +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; + + Error[n] = error; psi = 2.0*(f0*(1.0 - rlx) + rlx*(sum_q + 0.125*rho_e)); - - idx = Map[n]; - Psi[idx] = psi; - + + idx = Map[n]; + Psi[idx] = psi; + // q = 0 dist[n] = W0*psi;// @@ -593,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; - //........................................................................ } } @@ -635,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, 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)); + } +} +// +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, Np); + 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, 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)); + } +} +// +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, Np); + 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/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 diff --git a/hip/Poisson.hip b/hip/Poisson.hip index 2808ffa1..ac36000d 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>>(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, 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, 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, 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, diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 21f55225..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"); @@ -1549,7 +1550,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 4dab7323..ee4ccf1e 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" ); } @@ -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,18 +797,17 @@ 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(); //************************************************************************/ - // Check convergence of steady-state solution if (timestep==2){ //save electric potential for convergence check @@ -863,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,timestep);//perform collision + ScaLBL_Comm->Barrier(); comm.barrier(); + + // *************EVEN TIMESTEP*************// + timestep++; + //SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential + SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC,timestep);//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); @@ -1008,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); @@ -1029,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(); @@ -1038,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); @@ -1056,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(); @@ -1182,6 +1369,53 @@ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){ } } +void ScaLBL_Poisson::WriteVis( int timestep) { + + auto vis_db = db->getDatabase("Visualization"); + 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 8ab60140..e712f76e 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); @@ -35,12 +36,15 @@ 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, DoubleArray &Values_z); void getElectricField_debug(int timestep); void Checkpoint(); + void WriteVis( int timestep); void DummyChargeDensity(); //for debugging @@ -114,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/sample_scripts/configure_arden b/sample_scripts/configure_arden index 89ffb766..cdf524cf 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 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"); diff --git a/tests/lbpm_nernst_planck_cell_simulator.cpp b/tests/lbpm_nernst_planck_cell_simulator.cpp index 3ac2065c..d038b782 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 @@ -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); 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); //}