From e6d74f12d2f36b306e9cb360c89e88aaae10121e Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 9 Jan 2016 21:37:44 -0500 Subject: [PATCH] Updated TwoPhase.cpp with additional kinematic quantities for common curve dynamics --- common/TwoPhase.cpp | 29 +++++++++++++++++------------ common/TwoPhase.h | 2 ++ common/pmmc.h | 11 +++++++++-- 3 files changed, 28 insertions(+), 14 deletions(-) diff --git a/common/TwoPhase.cpp b/common/TwoPhase.cpp index da2b4320..845660ea 100644 --- a/common/TwoPhase.cpp +++ b/common/TwoPhase.cpp @@ -149,7 +149,9 @@ TwoPhase::TwoPhase(Domain &dm): fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz "); 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+"); @@ -450,8 +452,15 @@ void TwoPhase::ComputeLocal() // wn interface averages if (n_nw_pts > 0){ 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 0){ 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); - 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); - pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y, SDs_z, KNwns_values, KGwns_values, KNwns, KGwns, 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(&efawns,&efawns_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 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); @@ -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 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 - 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 fflush(TIMELOG); } diff --git a/common/TwoPhase.h b/common/TwoPhase.h index fc2f38b3..3b2e373f 100644 --- a/common/TwoPhase.h +++ b/common/TwoPhase.h @@ -98,6 +98,8 @@ public: double wp_volume_global; // volume for the wetting phase double As_global; double wwndnw, wwndnw_global; + double wwnsdnwn, wwnsdnwn_global; + double Jwnwwndnw, Jwnwwndnw_global; double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0) DoubleArray van; DoubleArray vaw; diff --git a/common/pmmc.h b/common/pmmc.h index 746d9a8e..87cd48db 100644 --- a/common/pmmc.h +++ b/common/pmmc.h @@ -3969,7 +3969,7 @@ inline double pmmc_CubeSurfaceOrientation(DoubleArray &Orientation, DTMutableLis 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 &Sx, DoubleArray &Sy, DoubleArray &Sz, DTMutableList &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; Point P,A,B; lwns = 0.0; + double ReturnValue = 0; NULL_USE(lwns); 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(1) += zeta*nwns_y*s; ReturnVector(2) += zeta*nwns_z*s; + ReturnValue += zeta*s; } } NULL_USE(tangent_x); NULL_USE(tangent_y); NULL_USE(tangent_z); + + return ReturnValue; } inline void pmmc_CurveOrientation(DoubleArray &Orientation, DTMutableList &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 &Points, IntArray &Triangles, DoubleArray &SurfaceVector, DoubleArray &AvgSpeed, DoubleArray &AvgVel, 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 s,s1,s2,s3,area; double norm, zeta; + double ReturnValue=0.0; TriLinPoly Px,Py,Pz,Pt; 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(2) += area*zeta*z; AvgSpeed(r) = zeta*area; + ReturnValue += zeta*area; } } + return ReturnValue; //............................................................................. } //--------------------------------------------------------------------------------------------------------