Skip to content

spurious diffusive effect correction #112

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Ricardoleite wants to merge 2 commits into
base: dev
Choose a base branch
from
+272 −33
Open

spurious diffusive effect correction #112

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
474d3c2
spurious diffusive effect correction
Ricardoleite Feb 20, 2026
656879c
missing ID_solid
Ricardoleite Feb 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions common/ScaLBL.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);

Expand Down
124 changes: 116 additions & 8 deletions cpu/Color.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -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;
Expand Down Expand Up @@ -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));
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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) {

Expand All @@ -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;
Expand Down Expand Up @@ -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));
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down
Loading
Loading