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

Updated TwoPhase.cpp with additional kinematic quantities for common curve dynamics

Browse Source
This commit is contained in:
James E McClure 2016-01-09 21:37:44 -05:00
parent 45886887bc
commit e6d74f12d2
3 changed files with 28 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
common/TwoPhase.cpp
View File

@ -149,7 +149,9 @@ TwoPhase::TwoPhase(Domain &dm):
fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz "); fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz "); fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
fprintf(TIMELOG,"trawn trJwn trRwn wwndnw Euler Kn Jn An\n"); // trimmed curvature & minkowski measures fprintf(TIMELOG,"trawn trJwn trRwn "); //trimmed curvature,
fprintf(TIMELOG,"wwndnw wwnsdnwn Jwnwwndnw "); //kinematic quantities,
fprintf(TIMELOG,"Euler Kn Jn An\n"); //miknowski measures,
} }
NWPLOG = fopen("components.NWP.tcat","a+"); NWPLOG = fopen("components.NWP.tcat","a+");
@ -450,8 +452,15 @@ void TwoPhase::ComputeLocal()
// wn interface averages // wn interface averages
if (n_nw_pts > 0){ if (n_nw_pts > 0){
awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris); awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris); Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
// Compute the normal speed of the interface
wwndnw += pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
//for (int p=0; p <n_nw_tris; p++) wwndnw += InterfaceSpeed(p);
for (int p=0; p <n_nw_tris; p++) Jwnwwndnw += InterfaceSpeed(p)*Values(p);
// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels) // Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
@ -460,22 +469,15 @@ void TwoPhase::ComputeLocal()
pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues, pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SDs,nw_pts,nw_tris,Values,DistanceValues,
i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn); i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn);
// Compute the normal speed of the interface
pmmc_InterfaceSpeed(dPdt, SDn_x, SDn_y, SDn_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
for (int p=0; p <n_nw_tris; p++) wwndnw += InterfaceSpeed(p);
} }
// wns common curve averages // wns common curve averages
if (n_local_nws_pts > 0){ if (n_local_nws_pts > 0){
efawns += pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z, efawns += pmmc_CubeContactAngle(CubeValues,Values,SDn_x,SDn_y,SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts); local_nws_pts,i,j,k,n_local_nws_pts);
pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z, wwnsdnwn += pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z,SDs_x,SDs_y,SDs_z,
local_nws_pts,i,j,k,n_local_nws_pts); local_nws_pts,i,j,k,n_local_nws_pts);
pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y, pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
SDs_z, KNwns_values, KGwns_values, KNwns, KGwns, SDs_z, KNwns_values, KGwns_values, KNwns, KGwns,
nws_pts, n_nws_pts, i, j, k); nws_pts, n_nws_pts, i, j, k);
@ -1107,6 +1109,8 @@ void TwoPhase::Reduce()
MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); MPI_Allreduce(&KNwns,&KNwns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wwndnw,&wwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); MPI_Allreduce(&wwndnw,&wwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&wwnsdnwn,&wwnsdnwn_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&Jwnwwndnw,&Jwnwwndnw_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
// Phase averages // Phase averages
MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm); MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
@ -1212,7 +1216,8 @@ void TwoPhase::PrintAll(int timestep)
Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ", fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global, wwndnw_global); // Trimmed curvature fprintf(TIMELOG,"%.5g %.5g %.5g ",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature
fprintf(TIMELOG,"%.5g %.5g %.5g ",wwndnw_global, wwnsdnwn_global, Jwnwwndnw_global); // kinematic quantities
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global); // minkowski measures fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global); // minkowski measures
fflush(TIMELOG); fflush(TIMELOG);
} }

2
common/TwoPhase.h
View File

@ -98,6 +98,8 @@ public:
double wp_volume_global; // volume for the wetting phase double wp_volume_global; // volume for the wetting phase
double As_global; double As_global;
double wwndnw, wwndnw_global; double wwndnw, wwndnw_global;
double wwnsdnwn, wwnsdnwn_global;
double Jwnwwndnw, Jwnwwndnw_global;
double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0) double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0)
DoubleArray van; DoubleArray van;
DoubleArray vaw; DoubleArray vaw;

11
common/pmmc.h
View File

@ -3969,7 +3969,7 @@ inline double pmmc_CubeSurfaceOrientation(DoubleArray &Orientation, DTMutableLis
return area; return area;
} }
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, DoubleArray &ReturnVector, inline double pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, DoubleArray &ReturnVector,
DoubleArray &Fx, DoubleArray &Fy, DoubleArray &Fz, DoubleArray &Fx, DoubleArray &Fy, DoubleArray &Fz,
DoubleArray &Sx, DoubleArray &Sy, DoubleArray &Sz, DoubleArray &Sx, DoubleArray &Sy, DoubleArray &Sz,
DTMutableList<Point> &Points, int i, int j, int k, int npts) DTMutableList<Point> &Points, int i, int j, int k, int npts)
@ -3983,6 +3983,7 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
double nwns_x, nwns_y, nwns_z; double nwns_x, nwns_y, nwns_z;
Point P,A,B; Point P,A,B;
lwns = 0.0; lwns = 0.0;
double ReturnValue = 0;
NULL_USE(lwns); NULL_USE(lwns);
TriLinPoly Px,Py,Pz,SDx,SDy,SDz,Pt; TriLinPoly Px,Py,Pz,SDx,SDy,SDz,Pt;
@ -4074,11 +4075,14 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
ReturnVector(0) += zeta*nwns_x*s; ReturnVector(0) += zeta*nwns_x*s;
ReturnVector(1) += zeta*nwns_y*s; ReturnVector(1) += zeta*nwns_y*s;
ReturnVector(2) += zeta*nwns_z*s; ReturnVector(2) += zeta*nwns_z*s;
ReturnValue += zeta*s;
} }
} }
NULL_USE(tangent_x); NULL_USE(tangent_x);
NULL_USE(tangent_y); NULL_USE(tangent_y);
NULL_USE(tangent_z); NULL_USE(tangent_z);
return ReturnValue;
} }
inline void pmmc_CurveOrientation(DoubleArray &Orientation, DTMutableList<Point> &Points, int npts, int i, int j, int k){ inline void pmmc_CurveOrientation(DoubleArray &Orientation, DTMutableList<Point> &Points, int npts, int i, int j, int k){
@ -4315,7 +4319,7 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------
inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray &P_y, DoubleArray &P_z, inline double pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray &P_y, DoubleArray &P_z,
DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles, DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &SurfaceVector, DoubleArray &AvgSpeed, DoubleArray &AvgVel, DoubleArray &SurfaceVector, DoubleArray &AvgSpeed, DoubleArray &AvgVel,
int i, int j, int k, int npts, int ntris) int i, int j, int k, int npts, int ntris)
@ -4324,6 +4328,7 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
double x,y,z; double x,y,z;
double s,s1,s2,s3,area; double s,s1,s2,s3,area;
double norm, zeta; double norm, zeta;
double ReturnValue=0.0;
TriLinPoly Px,Py,Pz,Pt; TriLinPoly Px,Py,Pz,Pt;
Px.assign(P_x,i,j,k); Px.assign(P_x,i,j,k);
@ -4365,8 +4370,10 @@ inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray
AvgVel(1) += area*zeta*y; AvgVel(1) += area*zeta*y;
AvgVel(2) += area*zeta*z; AvgVel(2) += area*zeta*z;
AvgSpeed(r) = zeta*area; AvgSpeed(r) = zeta*area;
ReturnValue += zeta*area;
} }
} }
return ReturnValue;
//............................................................................. //.............................................................................
} }
//-------------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------------