diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 9addb648..74ec2e0d 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -628,7 +628,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_MRT(int *neighborList, double *dist, * @param Np - size of local sub-domain (derived from Domain structure) */ extern "C" void ScaLBL_D3Q19_AAeven_Color( - int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, signed char *IDSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); @@ -660,7 +660,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( */ extern "C" void ScaLBL_D3Q19_AAodd_Color( int *NeighborList, int *Map, double *dist, double *Aq, double *Bq, - double *Den, double *Phi, double *Vel, double rhoA, double rhoB, + double *Den, double *Phi, signed char *IDSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); diff --git a/cpu/Color.cpp b/cpu/Color.cpp index 9e32147c..43d28c29 100644 --- a/cpu/Color.cpp +++ b/cpu/Color.cpp @@ -1429,7 +1429,7 @@ extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, // double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, // double Fx, double Fy, double Fz, int start, int finish, int Np){ extern "C" void ScaLBL_D3Q19_AAeven_Color( - int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, + int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, signed char *IDSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np) { @@ -1446,6 +1446,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( double C, nx, ny, nz; //color gradient magnitude and direction double ux, uy, uz; double phi, tau, rho0, rlx_setA, rlx_setB; + signed char id1, id2, id3, id4, id5, id6, id7, id8, id9, id10, id11, id12, id13, id14, id15, id16, id17, id18; + double nsx, nsy, nsz; // Rock Fluid interface normal vector + double npx, npy, npz; // contact angle vector const double mrt_V1 = 0.05263157894736842; const double mrt_V2 = 0.012531328320802; @@ -1484,57 +1487,75 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //........................................................................ nn = ijk - 1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 + id1 = IDSolid[nn]; //........................................................................ nn = ijk + 1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 + id2 = IDSolid[nn]; //........................................................................ nn = ijk - strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 + id3 = IDSolid[nn]; //........................................................................ nn = ijk + strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 + id4 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 + id5 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 + id6 = IDSolid[nn]; //........................................................................ nn = ijk - strideY - 1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 + id7 = IDSolid[nn]; //........................................................................ nn = ijk + strideY + 1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 + id8 = IDSolid[nn]; //........................................................................ nn = ijk + strideY - 1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 + id9 = IDSolid[nn]; //........................................................................ nn = ijk - strideY + 1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 + id10 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ - 1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 + id11 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ + 1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 + id12 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ - 1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 + id13 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ + 1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 + id14 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ - strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 + id15 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ + strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 + id16 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ - strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 + id17 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ + strideY; // neighbor index (get convention) m18 = Phi[nn]; // get neighbor for phi - 18 + id18 = IDSolid[nn]; //............Compute the Color Gradient................................... nx = -(m1 - m2 + 0.5 * (m7 - m8 + m9 - m10 + m11 - m12 + m13 - m14)); ny = -(m3 - m4 + 0.5 * (m7 - m8 - m9 + m10 + m15 - m16 + m17 - m18)); @@ -1549,6 +1570,39 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( ny = ny / ColorMag; nz = nz / ColorMag; + //...........Correct wettability vector for Mass Balance................................. + if (id1 == 0 || id2 == 0 || id3 == 0 || id4 == 0 || id5 == 0 || id6 == 0 || id7 == 0 || id8 == 0 || id9 == 0 || id10 == 0 || + id11 == 0 || id12 == 0 || id13 == 0 || id14 == 0 || id15 == 0 || id16 == 0 || id17 == 0 || id18 == 0){ + + int int_nsx = (id1 - id2) * 2 + (id7 - id8 + id9 - id10 + id11 - id12 + id13 - id14); + int int_nsy = (id3 - id4) * 2 + (id7 - id8 - id9 + id10 + id15 - id16 + id17 - id18); + int int_nsz = (id5 - id6) * 2 + (id11 - id12 - id13 + id14 + id15 - id16 - id17 + id18); + + double Mag = sqrt(double(int_nsx * int_nsx + int_nsy * int_nsy + int_nsz * int_nsz)); + if (Mag == 0.0) + Mag = 1.0f; + nsx = -int_nsx / Mag; + nsy = -int_nsy / Mag; + nsz = -int_nsz / Mag; + + int countid = id1 + id2 + id3 + id4 + id5 + id6 + id7 + id8 + id9 + id10 + id11 + id12 + id13 + id14 + id15 + id16 + id17 + id18; + + double aff = (m1*(id1-1) + m2*(id2-1) + m3*(id3-1) + m4*(id4-1) + m5*(id5-1) + m6*(id6-1) + m7*(id7-1) + m8*(id8-1) + + m9*(id9-1) + m10*(id10-1) + m11*(id11-1) + m12*(id12-1) + m13*(id13-1) + m14*(id14-1) + m15*(id15-1) + + m16*(id16-1) + m17*(id17-1) + m18*(id18-1)) / (18.0f - double(countid)); + + Mag = 1.0f-(nx*nsx + ny*nsy + nz*nsz)*(nx*nsx + ny*nsy + nz*nsz); + Mag = (Mag > 0.0f) ? sqrtf(Mag) : 1.0f; + + npx = (nx - nsx*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsx*aff; + npy = (ny - nsy*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsy*aff; + npz = (nz - nsz*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsz*aff; + }else{ + npx = nx; + npy = ny; + npz = nz; + } + // q=0 fq = dist[n]; rho = fq; @@ -1989,7 +2043,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nx; + delta = beta * nA * nB * nAB * 0.1111111111111111 * npx; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * ux)) + delta; @@ -2005,7 +2059,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 2 // Cq = {0,1,0} - delta = beta * nA * nB * nAB * 0.1111111111111111 * ny; + delta = beta * nA * nB * nAB * 0.1111111111111111 * npy; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uy)) + delta; @@ -2020,7 +2074,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( //............................................... // q = 4 // Cq = {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nz; + delta = beta * nA * nB * nAB * 0.1111111111111111 * npz; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uz)) + delta; @@ -2041,7 +2095,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color( // double Fx, double Fy, double Fz, int start, int finish, int Np){ extern "C" void ScaLBL_D3Q19_AAodd_Color( int *neighborList, int *Map, double *dist, double *Aq, double *Bq, - double *Den, double *Phi, double *Vel, double rhoA, double rhoB, + double *Den, double *Phi, signed char *IDSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np) { @@ -2061,6 +2115,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( double C, nx, ny, nz; //color gradient magnitude and direction double ux, uy, uz; double phi, tau, rho0, rlx_setA, rlx_setB; + signed char id1, id2, id3, id4, id5, id6, id7, id8, id9, id10, id11, id12, id13, id14, id15, id16, id17, id18; + double nsx, nsy, nsz; // Rock Fluid interface normal vector + double npx, npy, npz; // contact angle vector const double mrt_V1 = 0.05263157894736842; const double mrt_V2 = 0.012531328320802; @@ -2099,57 +2156,75 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //........................................................................ nn = ijk - 1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 + id1 = IDSolid[nn]; // 0 if Rock 1 if Fluid //........................................................................ nn = ijk + 1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 + id2 = IDSolid[nn]; //........................................................................ nn = ijk - strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 + id3 = IDSolid[nn]; //........................................................................ nn = ijk + strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 + id4 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 + id5 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 + id6 = IDSolid[nn]; //........................................................................ nn = ijk - strideY - 1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 + id7 = IDSolid[nn]; //........................................................................ nn = ijk + strideY + 1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 + id8 = IDSolid[nn]; //........................................................................ nn = ijk + strideY - 1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 + id9 = IDSolid[nn]; //........................................................................ nn = ijk - strideY + 1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 + id10 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ - 1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 + id11 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ + 1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 + id12 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ - 1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 + id13 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ + 1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 + id14 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ - strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 + id15 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ + strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 + id16 = IDSolid[nn]; //........................................................................ nn = ijk + strideZ - strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 + id17 = IDSolid[nn]; //........................................................................ nn = ijk - strideZ + strideY; // neighbor index (get convention) m18 = Phi[nn]; // get neighbor for phi - 18 + id18 = IDSolid[nn]; //............Compute the Color Gradient................................... nx = -(m1 - m2 + 0.5 * (m7 - m8 + m9 - m10 + m11 - m12 + m13 - m14)); ny = -(m3 - m4 + 0.5 * (m7 - m8 - m9 + m10 + m15 - m16 + m17 - m18)); @@ -2164,6 +2239,39 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( ny = ny / ColorMag; nz = nz / ColorMag; + //...........Correct wettability vector................................. + if (id1 == 0 || id2 == 0 || id3 == 0 || id4 == 0 || id5 == 0 || id6 == 0 || id7 == 0 || id8 == 0 || id9 == 0 || id10 == 0 || + id11 == 0 || id12 == 0 || id13 == 0 || id14 == 0 || id15 == 0 || id16 == 0 || id17 == 0 || id18 == 0){ + + int int_nsx = (id1 - id2) * 2 + (id7 - id8 + id9 - id10 + id11 - id12 + id13 - id14); + int int_nsy = (id3 - id4) * 2 + (id7 - id8 - id9 + id10 + id15 - id16 + id17 - id18); + int int_nsz = (id5 - id6) * 2 + (id11 - id12 - id13 + id14 + id15 - id16 - id17 + id18); + + double Mag = sqrt(double(int_nsx * int_nsx + int_nsy * int_nsy + int_nsz * int_nsz)); + if (Mag == 0.0) + Mag = 1.0f; + nsx = -int_nsx / Mag; + nsy = -int_nsy / Mag; + nsz = -int_nsz / Mag; + + int countid = id1 + id2 + id3 + id4 + id5 + id6 + id7 + id8 + id9 + id10 + id11 + id12 + id13 + id14 + id15 + id16 + id17 + id18; + + double aff = (m1*(id1-1) + m2*(id2-1) + m3*(id3-1) + m4*(id4-1) + m5*(id5-1) + m6*(id6-1) + m7*(id7-1) + m8*(id8-1) + + m9*(id9-1) + m10*(id10-1) + m11*(id11-1) + m12*(id12-1) + m13*(id13-1) + m14*(id14-1) + m15*(id15-1) + + m16*(id16-1) + m17*(id17-1) + m18*(id18-1)) / (18.0f - double(countid)); + + Mag = 1.0f-(nx*nsx + ny*nsy + nz*nsz)*(nx*nsx + ny*nsy + nz*nsz); + Mag = (Mag > 0.0f) ? sqrtf(Mag) : 1.0f; + + npx = (nx - nsx*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsx*aff; + npy = (ny - nsy*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsy*aff; + npz = (nz - nsz*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsz*aff; + }else{ + npx = nx; + npy = ny; + npz = nz; + } + // q=0 fq = dist[n]; rho = fq; @@ -2666,7 +2774,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nx; + delta = beta * nA * nB * nAB * 0.1111111111111111 * npx; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * ux)) + delta; @@ -2685,7 +2793,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // Cq = {0,1,0} - delta = beta * nA * nB * nAB * 0.1111111111111111 * ny; + delta = beta * nA * nB * nAB * 0.1111111111111111 * npy; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uy)) + delta; @@ -2705,7 +2813,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color( //............................................... // q = 4 // Cq = {0,0,1} - delta = beta * nA * nB * nAB * 0.1111111111111111 * nz; + delta = beta * nA * nB * nAB * 0.1111111111111111 * npz; if (!(nA * nB * nAB > 0)) delta = 0; a1 = nA * (0.1111111111111111 * (1 + 4.5 * uz)) + delta; diff --git a/cuda/Color.cu b/cuda/Color.cu index fc35c5c9..23ad547e 100644 --- a/cuda/Color.cu +++ b/cuda/Color.cu @@ -1266,7 +1266,7 @@ __global__ void dvc_ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int } -__global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, +__global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, signed char *IDSolid, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; @@ -1281,6 +1281,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; + signed char id1, id2, id3, id4, id5, id6, id7, id8, id9, id10, id11, id12, id13, id14, id15, id16, id17, id18; + double nsx, nsy, nsz; // Rock Fluid interface normal vector + double npx, npy, npz; // contact angle vector const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1323,57 +1326,75 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //........................................................................ nn = ijk-1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 + id1 = IDSolid[nn]; //........................................................................ nn = ijk+1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 + id2 = IDSolid[nn]; //........................................................................ nn = ijk-strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 + id3 = IDSolid[nn]; //........................................................................ nn = ijk+strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 + id4 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 + id5 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 + id6 = IDSolid[nn]; //........................................................................ nn = ijk-strideY-1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 + id7 = IDSolid[nn]; //........................................................................ nn = ijk+strideY+1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 + id8 = IDSolid[nn]; //........................................................................ nn = ijk+strideY-1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 + id9 = IDSolid[nn]; //........................................................................ nn = ijk-strideY+1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 + id10 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ-1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 + id11 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ+1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 + id12 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ-1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 + id13 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ+1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 + id14 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ-strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 + id15 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ+strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 + id16 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ-strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 + id17 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ+strideY; // neighbor index (get convention) m18 = Phi[nn]; // get neighbor for phi - 18 + id18 = IDSolid[nn]; //............Compute the Color Gradient................................... nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); @@ -1387,6 +1408,40 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A ny = ny/ColorMag; nz = nz/ColorMag; + //...........Correct wettability vector for Mass Balance................................. + if (id1 == 0 || id2 == 0 || id3 == 0 || id4 == 0 || id5 == 0 || id6 == 0 || id7 == 0 || id8 == 0 || id9 == 0 || id10 == 0 || + id11 == 0 || id12 == 0 || id13 == 0 || id14 == 0 || id15 == 0 || id16 == 0 || id17 == 0 || id18 == 0){ + + int int_nsx = (id1 - id2) * 2 + (id7 - id8 + id9 - id10 + id11 - id12 + id13 - id14); + int int_nsy = (id3 - id4) * 2 + (id7 - id8 - id9 + id10 + id15 - id16 + id17 - id18); + int int_nsz = (id5 - id6) * 2 + (id11 - id12 - id13 + id14 + id15 - id16 - id17 + id18); + + double Mag = sqrt(double(int_nsx * int_nsx + int_nsy * int_nsy + int_nsz * int_nsz)); + if (Mag == 0.0) + Mag = 1.0f; + nsx = -int_nsx / Mag; + nsy = -int_nsy / Mag; + nsz = -int_nsz / Mag; + + int countid = id1 + id2 + id3 + id4 + id5 + id6 + id7 + id8 + id9 + id10 + id11 + id12 + id13 + id14 + id15 + id16 + id17 + id18; + + double aff = (m1*(id1-1) + m2*(id2-1) + m3*(id3-1) + m4*(id4-1) + m5*(id5-1) + m6*(id6-1) + m7*(id7-1) + m8*(id8-1) + + m9*(id9-1) + m10*(id10-1) + m11*(id11-1) + m12*(id12-1) + m13*(id13-1) + m14*(id14-1) + m15*(id15-1) + + m16*(id16-1) + m17*(id17-1) + m18*(id18-1)) / (18.0f - double(countid)); + + Mag = 1.0f-(nx*nsx + ny*nsy + nz*nsz)*(nx*nsx + ny*nsy + nz*nsz); + Mag = (Mag > 0.0f) ? sqrtf(Mag) : 1.0f; + + npx = (nx - nsx*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsx*aff; + npy = (ny - nsy*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsy*aff; + npz = (nz - nsz*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsz*aff; + }else{ + npx = nx; + npy = ny; + npz = nz; + } + + // q=0 fq = dist[n]; rho = fq; @@ -1802,7 +1857,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nx; + delta = beta*nA*nB*nAB*0.1111111111111111*npx; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -1817,7 +1872,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 2 // Cq = {0,1,0} - delta = beta*nA*nB*nAB*0.1111111111111111*ny; + delta = beta*nA*nB*nAB*0.1111111111111111*npy; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -1831,7 +1886,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A //............................................... // q = 4 // Cq = {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nz; + delta = beta*nA*nB*nAB*0.1111111111111111*npz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -1850,7 +1905,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double *Phi, signed char *IDSolid, double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ int n,nn,ijk,nread; @@ -1869,6 +1924,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double double C,nx,ny,nz; //color gradient magnitude and direction double ux,uy,uz; double phi,tau,rho0,rlx_setA,rlx_setB; + signed char id1, id2, id3, id4, id5, id6, id7, id8, id9, id10, id11, id12, id13, id14, id15, id16, id17, id18; + double nsx, nsy, nsz; // Rock Fluid interface normal vector + double npx, npy, npz; // contact angle vector const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1910,57 +1968,75 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //........................................................................ nn = ijk-1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 + id1 = IDSolid[nn]; // 0 if Rock 1 if Fluid //........................................................................ nn = ijk+1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 + id2 = IDSolid[nn]; //........................................................................ nn = ijk-strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 + id3 = IDSolid[nn]; //........................................................................ nn = ijk+strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 + id4 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 + id5 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 + id6 = IDSolid[nn]; //........................................................................ nn = ijk-strideY-1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 + id7 = IDSolid[nn]; //........................................................................ nn = ijk+strideY+1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 + id8 = IDSolid[nn]; //........................................................................ nn = ijk+strideY-1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 + id9 = IDSolid[nn]; //........................................................................ nn = ijk-strideY+1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 + id10 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ-1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 + id11 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ+1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 + id12 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ-1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 + id13 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ+1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 + id14 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ-strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 + id15 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ+strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 + id16 = IDSolid[nn]; //........................................................................ nn = ijk+strideZ-strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 + id17 = IDSolid[nn]; //........................................................................ nn = ijk-strideZ+strideY; // neighbor index (get convention) m18 = Phi[nn]; // get neighbor for phi - 18 + id18 = IDSolid[nn]; //............Compute the Color Gradient................................... nx = -(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); ny = -(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); @@ -1972,7 +2048,40 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double if (C==0.0) ColorMag=1.0; nx = nx/ColorMag; ny = ny/ColorMag; - nz = nz/ColorMag; + nz = nz/ColorMag; + + //...........Correct wettability vector for Mass Balance................................. + if (id1 == 0 || id2 == 0 || id3 == 0 || id4 == 0 || id5 == 0 || id6 == 0 || id7 == 0 || id8 == 0 || id9 == 0 || id10 == 0 || + id11 == 0 || id12 == 0 || id13 == 0 || id14 == 0 || id15 == 0 || id16 == 0 || id17 == 0 || id18 == 0){ + + int int_nsx = (id1 - id2) * 2 + (id7 - id8 + id9 - id10 + id11 - id12 + id13 - id14); + int int_nsy = (id3 - id4) * 2 + (id7 - id8 - id9 + id10 + id15 - id16 + id17 - id18); + int int_nsz = (id5 - id6) * 2 + (id11 - id12 - id13 + id14 + id15 - id16 - id17 + id18); + + double Mag = sqrt(double(int_nsx * int_nsx + int_nsy * int_nsy + int_nsz * int_nsz)); + if (Mag == 0.0) + Mag = 1.0f; + nsx = -int_nsx / Mag; + nsy = -int_nsy / Mag; + nsz = -int_nsz / Mag; + + int countid = id1 + id2 + id3 + id4 + id5 + id6 + id7 + id8 + id9 + id10 + id11 + id12 + id13 + id14 + id15 + id16 + id17 + id18; + + double aff = (m1*(id1-1) + m2*(id2-1) + m3*(id3-1) + m4*(id4-1) + m5*(id5-1) + m6*(id6-1) + m7*(id7-1) + m8*(id8-1) + + m9*(id9-1) + m10*(id10-1) + m11*(id11-1) + m12*(id12-1) + m13*(id13-1) + m14*(id14-1) + m15*(id15-1) + + m16*(id16-1) + m17*(id17-1) + m18*(id18-1)) / (18.0f - double(countid)); + + Mag = 1.0f-(nx*nsx + ny*nsy + nz*nsz)*(nx*nsx + ny*nsy + nz*nsz); + Mag = (Mag > 0.0f) ? sqrtf(Mag) : 1.0f; + + npx = (nx - nsx*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsx*aff; + npy = (ny - nsy*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsy*aff; + npz = (nz - nsz*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsz*aff; + }else{ + npx = nx; + npy = ny; + npz = nz; + } // q=0 fq = dist[n]; @@ -2451,7 +2560,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nx; + delta = beta*nA*nB*nAB*0.1111111111111111*npx; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2469,7 +2578,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // Cq = {0,1,0} - delta = beta*nA*nB*nAB*0.1111111111111111*ny; + delta = beta*nA*nB*nAB*0.1111111111111111*npy; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2488,7 +2597,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double //............................................... // q = 4 // Cq = {0,0,1} - delta = beta*nA*nB*nAB*0.1111111111111111*nz; + delta = beta*nA*nB*nAB*0.1111111111111111*npz; if (!(nA*nB*nAB>0)) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -4007,14 +4116,14 @@ extern "C" void ScaLBL_D3Q7_ColorCollideMass(char *ID, double *A_even, double *A } // Pressure Boundary Conditions Functions -extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, +extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi, signed char *IDSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ cudaProfilerStart(); cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_Color, cudaFuncCachePreferL1); - dvc_ScaLBL_D3Q19_AAeven_Color<<>>(Map, dist, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB, + dvc_ScaLBL_D3Q19_AAeven_Color<<>>(Map, dist, Aq, Bq, Den, Phi, IDSolid, Vel, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -4025,13 +4134,13 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do } extern "C" void ScaLBL_D3Q19_AAodd_Color(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, + double *Phi, signed char *IDSolid, double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ cudaProfilerStart(); cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_Color, cudaFuncCachePreferL1); - dvc_ScaLBL_D3Q19_AAodd_Color<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, Vel, + dvc_ScaLBL_D3Q19_AAodd_Color<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, IDSolid, Vel, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 88caa2ea..64899510 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -32,7 +32,7 @@ ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, tauB(0), rhoA(0), rhoB(0), alpha(0), beta(0), Fx(0), Fy(0), Fz(0), flux(0), din(0), dout(0), inletA(0), inletB(0), outletA(0), outletB(0), Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), - BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), + BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), IDSolid(nullptr), NeighborList(nullptr), dvcMap(nullptr), fq(nullptr), Aq(nullptr), Bq(nullptr), Den(nullptr), Phi(nullptr), ColorGrad(nullptr), Velocity(nullptr), Pressure(nullptr), comm(COMM) { @@ -448,6 +448,7 @@ void ScaLBL_ColorModel::Create() { ScaLBL_AllocateDeviceMemory((void **)&Pressure, sizeof(double) * Np); ScaLBL_AllocateDeviceMemory((void **)&Velocity, 3 * sizeof(double) * Np); ScaLBL_AllocateDeviceMemory((void **)&ColorGrad, 3 * sizeof(double) * Np); + ScaLBL_AllocateDeviceMemory((void **)&IDSolid, sizeof(signed char) * Nx * Ny * Nz); //........................................................................... // Update GPU data structures if (rank == 0) @@ -503,6 +504,27 @@ void ScaLBL_ColorModel::Create() { if (rank == 0) printf("Model created \n"); delete[] PhaseLabel; + + // Initialize solid-fluid id + signed char *TmpID; + TmpID = new signed char[Nx * Ny * Nz]; + for (int k = 0; k < Nz; k++) { + for (int j = 0; j < Ny; j++) { + for (int i = 0; i < Nx; i++) { + int idx = k * Nx * Ny + j * Nx + i; + if (id[idx] <= 0){ + TmpID[idx] = 0; + } + else{ + TmpID[idx] = 1; + } + } + } + } + + ScaLBL_CopyToDevice(IDSolid, TmpID, sizeof(signed char) * Nx * Ny * Nz); + ScaLBL_Comm->Barrier(); + delete[] TmpID; } /******************************************************** @@ -690,7 +712,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); @@ -709,7 +731,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); @@ -735,7 +757,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -754,7 +776,7 @@ double ScaLBL_ColorModel::Run(int returntime) { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); @@ -1142,7 +1164,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm_Regular->SendHalo(Phi); ScaLBL_D3Q19_AAodd_Color( - NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm_Regular->RecvHalo(Phi); @@ -1161,7 +1183,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); @@ -1187,7 +1209,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); } ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); @@ -1206,7 +1228,7 @@ void ScaLBL_ColorModel::Run() { ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, IDSolid, Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); diff --git a/models/ColorModel.h b/models/ColorModel.h index 888f4fc2..625b6543 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -137,7 +137,7 @@ class ScaLBL_ColorModel { std::shared_ptr vis_db; IntArray Map; - signed char *id; + signed char *id, *IDSolid; int *NeighborList; int *dvcMap; double *fq, *Aq, *Bq;