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

intermediate contact angle nan fixes

Browse Source
This commit is contained in:
James McClure 2014-03-24 17:06:48 -04:00
parent 9eb4218a4e
commit 54936e6f68
3 changed files with 25 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
common/pmmc.h
View File

@ -3431,9 +3431,9 @@ inline int pmmc_CubeListFromMesh(IntArray &cubeList, int ncubes, int Nx, int Ny,
nc=0;
//...........................................................................
// Set up the cube list (very regular in this case due to lack of blob-ID)
for (k=0; k<Nz-2; k++){
for (j=0; j<Ny-2; j++){
for (i=0; i<Nx-2; i++){
for (k=1; k<Nz-2; k++){
for (j=1; j<Ny-2; j++){
for (i=1; i<Nx-2; i++){
cubeList(0,nc) = i;
cubeList(1,nc) = j;
cubeList(2,nc) = k;
@ -3645,13 +3645,14 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
// printf("grad F = %f \n", sqrt(pow(Fx(i,j,k),2)+pow(Fy(i,j,k),2)+pow(Fz(i,j,k),2)));
}
integral = 0.0;
@ -3664,8 +3665,12 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
// Compute the length of the segment
length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
integral += 0.5*length*(vA + vB);
/// printf("cos phi = %f, ", vA);
// printf("length = %f \n", length);
// printf("i,j,k = %i, %i, %i \n", i,j,k);
}
if (npts > 0)
// printf("integral =%f \n",integral);
return integral;
}
//--------------------------------------------------------------------------------------------------------

19
pmmc/Analysis.cpp
View File

@ -218,6 +218,7 @@ int main(int argc, char **argv)
fclose(VEL);
*/ //...........................................................................
// Calculate the time derivative of the phase indicator field
for (int n=0; n<Nx*Ny*Nz; n++) dPdt(n) = 0.5*(Phase_tplus(n) - Phase_tminus(n));
pmmc_MeshGradient(Phase,Phase_x,Phase_y,Phase_z,Nx,Ny,Nz);
@ -228,17 +229,16 @@ int main(int argc, char **argv)
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
if ( SignDist(i,j,k) > 0 ){
if ( SignDist(i,j,k) > 0.0 ){
// 1-D index for this cube corner
n = i + j*Nx + k*Nx*Ny;
// Compute the non-wetting phase volume contribution
if ( Phase(i,j,k) > 0.0 )
nwp_volume += 1.0;
// volume averages over the non-wetting phase
if ( Phase(i,j,k) > 0.9999 ){
if ( Phase(i,j,k) > 0.99 ){
// volume the excludes the interfacial region
vol_n += 1.0;
// pressure
@ -250,7 +250,7 @@ int main(int argc, char **argv)
}
// volume averages over the wetting phase
if ( Phase(i,j,k) < -0.9999 ){
if ( Phase(i,j,k) < -0.99 ){
// volume the excludes the interfacial region
vol_w += 1.0;
// pressure
@ -265,6 +265,9 @@ int main(int argc, char **argv)
}
}
printf("vol_n = %f \n",vol_n);
printf("vol_w = %f \n",vol_w);
// End of the loop to set the values
awn = aws = ans = lwns = 0.0;
nwp_volume = 0.0;
@ -290,6 +293,8 @@ int main(int argc, char **argv)
FILE *WNS_PTS;
WNS_PTS = fopen("wns-pts.out","w");
printf("ncubes = %i,\n",ncubes);
for (c=0;c<ncubes;c++){
// Get cube from the list
i = cubeList(0,c);
@ -322,7 +327,11 @@ int main(int argc, char **argv)
// Compute the average contact angle
efawns += pmmc_CubeContactAngle(CubeValues,ContactAngle,Phase_x,Phase_y,Phase_z,Sx,Sy,Sz,
local_nws_pts,i,j,k,n_local_nws_pts);
// printf("efawns= %f \n", efawns);
if (isnan(efawns))
c = ncubes;
// Compute the curvature of the wn interface
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature, nw_pts, nw_tris,
Curvature, i, j, k, n_nw_pts, n_nw_tris);

1
tests/lb2_Color_wia_mpi.cpp
View File

@ -1882,6 +1882,7 @@ int main(int argc, char **argv)
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
// Compute the gradients of the phase indicator and signed distance fields
pmmc_MeshGradient(Phase,Phase_x,Phase_y,Phase_z,Nx,Ny,Nz);
// pmmc_MeshGradient(SignDist,SignDist_x,SignDist_y,SignDist_z,Nx,Ny,Nz);