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

fix sign on geodesic curvature

Browse Source
This commit is contained in:
James McClure 2021-12-15 07:35:10 -05:00
parent d259434e5f
commit 5c794e0bd0
4 changed files with 53 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
analysis/TwoPhase.cpp
View File

@ -361,6 +361,7 @@ void TwoPhase::Initialize() {
wwndnw = 0.0;
wwnsdnwn = 0.0;
Jwnwwndnw = 0.0;
Xwn = Xns = Xws = 0.0;
}
void TwoPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB,
@ -780,6 +781,9 @@ void TwoPhase::ComputeStatic() {
Jwn += pmmc_CubeSurfaceInterpValue(
CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i,
j, k, n_nw_pts, n_nw_tris);
Xwn += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
n_nw_tris, i, j, k);
// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
pmmc_CubeTrimSurfaceInterpValues(
@ -812,6 +816,14 @@ void TwoPhase::ComputeStatic() {
lwns +=
pmmc_CubeCurveLength(local_nws_pts, n_local_nws_pts);
/* half contribution for vertices / edges at the common line
* each cube with contact line has a net of undercounting vertices
* each cube is undercounting edges due to internal counts
*/
Xwn += 0.25*n_local_nws_pts - 0.5;
Xws += 0.25*n_local_nws_pts - 0.5;
Xns += 0.25*n_local_nws_pts - 0.5;
}
// Solid interface averagees
@ -824,6 +836,12 @@ void TwoPhase::ComputeStatic() {
n_ns_tris);
aws += pmmc_CubeSurfaceOrientation(Gws, ws_pts, ws_tris,
n_ws_tris);
Xws += geomavg_EulerCharacteristic(ws_pts, ws_tris, n_ws_pts,
n_ws_tris, i, j, k);
Xns += geomavg_EulerCharacteristic(ns_pts, ns_tris, n_ns_pts,
n_ns_tris, i, j, k);
}
//...........................................................................
// Compute the integral curvature of the non-wetting phase
@ -848,9 +866,11 @@ void TwoPhase::ComputeStatic() {
Kn += pmmc_CubeSurfaceInterpValue(CubeValues, GaussCurvature,
nw_pts, nw_tris, Values, i, j,
k, n_nw_pts, n_nw_tris);
euler += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
n_nw_tris, i, j, k);
}
}
}
@ -1447,10 +1467,14 @@ void TwoPhase::Reduce() {
trawn_global = Dm->Comm.sumReduce(trawn);
trJwn_global = Dm->Comm.sumReduce(trJwn);
trRwn_global = Dm->Comm.sumReduce(trRwn);
euler_global = Dm->Comm.sumReduce(euler);
Xwn_global = Dm->Comm.sumReduce(Xwn);
Xws_global = Dm->Comm.sumReduce(Xws);
Xns_global = Dm->Comm.sumReduce(Xns);
An_global = Dm->Comm.sumReduce(An);
Jn_global = Dm->Comm.sumReduce(Jn);
Kn_global = Dm->Comm.sumReduce(Kn);
euler_global = Dm->Comm.sumReduce(euler);
Dm->Comm.barrier();
// Normalize the phase averages
@ -1533,7 +1557,7 @@ void TwoPhase::PrintStatic() {
if (fseek(STATIC, 0, SEEK_SET) == fseek(STATIC, 0, SEEK_CUR)) {
// If timelog is empty, write a short header to list the averages
fprintf(STATIC, "sw awn ans aws Jwn Kwn lwns cwns KGws "
"KGwn "); // Scalar averages
"KGwn Xwn Xws Xns "); // Scalar averages
fprintf(STATIC,
"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
fprintf(STATIC, "Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
@ -1553,6 +1577,8 @@ void TwoPhase::PrintStatic() {
fprintf(STATIC, "%.5g ", efawns_global); // average contact angle
fprintf(STATIC, "%.5g %.5g ", KNwns_global,
KGwns_global); // total curvature contributions of common line
fprintf(STATIC, "%.5g %.5g %.5g ", Xwn_global, Xns_global,
Xws_global); // Euler characteristic
fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gwn_global(0),
Gwn_global(1), Gwn_global(2), Gwn_global(3), Gwn_global(4),
Gwn_global(5)); // orientation of wn interface

2
analysis/TwoPhase.h
View File

@ -103,7 +103,9 @@ public:
double lwns_global;
double efawns, efawns_global; // averaged contact angle
double euler, Kn, Jn, An;
double Xwn, Xns, Xws;
double euler_global, Kn_global, Jn_global, An_global;
double Xwn_global, Xns_global, Xws_global;
double rho_n, rho_w;
double nu_n, nu_w;

20
analysis/pmmc.h
View File

@ -4572,14 +4572,13 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
nwsx /= norm;
nwsy /= norm;
nwsz /= norm;
/* JM does not remember why we should do this...
* if (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz < 0.0) {
nwnsx = -nwnsx;
nwnsy = -nwnsy;
nwnsz = -nwnsz;
/* normal to ws interface boundary should point into fluid (same direction as gradient) */
if (nwx * nwsx + nwy * nwsy + nwz * nwsz < 0.0) {
nwsx = -nwsx;
nwsy = -nwsy;
nwsz = -nwsz;
}
*/
// common curve normal in the fluid surface tangent plane (rel. geodesic curvature)
nwnx = twnsy * nwz - twnsz * nwy;
nwny = twnsz * nwx - twnsx * nwz;
@ -4590,6 +4589,13 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
nwnx /= norm;
nwny /= norm;
nwnz /= norm;
/* normal to wn interface boundary should point into the solid */
if (nsx * nwnx + nsy * nwny + nsz * nwnz > 0.0) {
nwnx = -nwnx;
nwny = -nwny;
nwnz = -nwnz;
}
if (length > 0.0) {
// normal curvature component in the direction of the solid surface

10
tests/lbpm_TwoPhase_analysis.cpp
View File

@ -126,6 +126,16 @@ int main(int argc, char **argv)
fillDouble.fill(Averages->SDs);
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(Averages->SDs,id,*Dm);
/* move the solid surface by a voxel to improve contact angle measurement
for (k=0;k<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
n=k*nx*ny+j*nx+i;
Averages->SDs(i,j,k) -= 1.0;
}
}
}*/
/* * * * * * * * * * * * * * * * * * * * * */
// Solve for the position of the fluid phase
for (k=0;k<nz;k++){