diff --git a/gpu/lb2_Color_wia_mpi.cpp b/gpu/lb2_Color_wia_mpi.cpp index a2d68614..45eb1f0a 100644 --- a/gpu/lb2_Color_wia_mpi.cpp +++ b/gpu/lb2_Color_wia_mpi.cpp @@ -9,7 +9,6 @@ using namespace std; - //************************************************************************* // Functions defined in Color.cu //************************************************************************* @@ -96,6 +95,7 @@ inline void PackID(int *list, int count, char *sendbuf, char *ID){ } } //*************************************************************************************** + inline void UnpackID(int *list, int count, char *recvbuf, char *ID){ // Fill in the phase ID values from neighboring processors // This unpacks the values once they have been recieved from neighbors @@ -106,7 +106,110 @@ inline void UnpackID(int *list, int count, char *recvbuf, char *ID){ ID[n] = recvbuf[idx]; } } + //*************************************************************************************** +inline void PackMeshData(int *list, int count, double *sendbuf, DoubleArray &Values){ + // Fill in the phase ID values from neighboring processors + // This packs up the values that need to be sent from one processor to another + int idx,n; + for (idx=0; idx tmp(20); DoubleArray Values(20); + DoubleArray InterfaceSpeed(20); + DoubleArray NormalVector(60); + DoubleArray vawn(3); // IntArray store; int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0; @@ -1922,6 +2063,7 @@ int main(int argc, char **argv) // Calculate the time derivative of the phase indicator field for (n=0; n 0.0){ temp = sqrt(temp); area += temp; - Orientation(0) += temp*nx*nx; // Gxx - Orientation(1) += temp*ny*ny; // Gyy - Orientation(2) += temp*nz*nz; // Gzz - Orientation(3) += temp*nx*ny; // Gxy - Orientation(4) += temp*nx*nz; // Gxz - Orientation(5) += temp*ny*nz; // Gyz + Orientation(0) += temp*nx*nx*normsq; // Gxx + Orientation(1) += temp*ny*ny*normsq; // Gyy + Orientation(2) += temp*nz*nz*normsq; // Gzz + Orientation(3) += temp*nx*ny*normsq; // Gxy + Orientation(4) += temp*nx*nz*normsq; // Gxz + Orientation(5) += temp*ny*nz*normsq; // Gyz } } return area; diff --git a/pmmc/Analysis.cpp b/pmmc/Analysis.cpp new file mode 100644 index 00000000..17e6d636 --- /dev/null +++ b/pmmc/Analysis.cpp @@ -0,0 +1,337 @@ +#include +#include +#include "pmmc.h" +#include "Domain.h" +//#include "PointList.h" +//#include "Array.h" + +#define CAPRAD 25 +#define RADIUS 20 +#define SPEED 1.0 + +using namespace std; + +int main(int argc, char **argv) +{ + //....................................................................... + // printf("Radius = %s \n,"RADIUS); + int Nx,Ny,Nz; + int i,j,k,p,q,r,n; + int nspheres; + double Lx,Ly,Lz; + //....................................................................... + Nx = Ny = Nz = 60; + //....................................................................... + // Reading the domain information file +/* //....................................................................... + ifstream domain("Domain.in"); + domain >> Nx; + domain >> Ny; + domain >> Nz; + domain >> nspheres; + domain >> Lx; + domain >> Ly; + domain >> Lz; +*/ //....................................................................... + + //....................................................................... + DoubleArray SignDist(Nx,Ny,Nz); + DoubleArray Phase(Nx,Ny,Nz); + DoubleArray Phase_x(Nx,Ny,Nz); + DoubleArray Phase_y(Nx,Ny,Nz); + DoubleArray Phase_z(Nx,Ny,Nz); + DoubleArray Sx(Nx,Ny,Nz); + DoubleArray Sy(Nx,Ny,Nz); + DoubleArray Sz(Nx,Ny,Nz); + DoubleArray Vel_x(Nx,Ny,Nz); + DoubleArray Vel_y(Nx,Ny,Nz); + DoubleArray Vel_z(Nx,Ny,Nz); + DoubleArray Press(Nx,Ny,Nz); + DoubleArray GaussCurvature(Nx,Ny,Nz); + DoubleArray MeanCurvature(Nx,Ny,Nz); + //....................................................................... + + //....................................................................... + double fluid_isovalue = 0.0; + double solid_isovalue = 0.0; + //....................................................................... + +/* //....................................................................... + double *cx,*cy,*cz,*rad; + cx = new double[nspheres]; + cy = new double[nspheres]; + cz = new double[nspheres]; + rad = new double[nspheres]; + //............................... + printf("Reading the sphere packing \n"); + ReadSpherePacking(nspheres,cx,cy,cz,rad); + //....................................................................... + //....................................................................... + // Compute the signed distance function for the sphere packing + SignedDistance(SignDist.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1); +*/ //....................................................................... + + /* **************************************************************** + VARIABLES FOR THE PMMC ALGORITHM + ****************************************************************** */ + //........................................................................... + // Averaging variables + //........................................................................... + double awn,ans,aws,lwns,nwp_volume; + double sw,vol_n,vol_w,paw,pan; + double efawns,Jwn; + double As; + double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0) + double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global; + double As_global; + // bool add=1; // Set to false if any corners contain nw-phase ( F > fluid_isovalue) + + int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0; + int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0; + + double s,s1,s2,s3; // Triangle sides (lengths) + Point A,B,C,P; + // double area; + int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}}; // cube corners + // int count_in=0,count_out=0; + // int nodx,nody,nodz; + // initialize lists for vertices for surfaces, common line + DTMutableList nw_pts(20); + DTMutableList ns_pts(20); + DTMutableList ws_pts(20); + DTMutableList nws_pts(20); + // initialize triangle lists for surfaces + IntArray nw_tris(3,20); + IntArray ns_tris(3,20); + IntArray ws_tris(3,20); + // initialize list for line segments + IntArray nws_seg(2,20); + DTMutableList tmp(20); + + // Initialize arrays for local solid surface + DTMutableList local_sol_pts(20); + int n_local_sol_pts = 0; + IntArray local_sol_tris(3,18); + int n_local_sol_tris; + DoubleArray values(20); + DTMutableList local_nws_pts(20); + int n_local_nws_pts; + + DoubleArray CubeValues(2,2,2); + DoubleArray ContactAngle(20); + DoubleArray Curvature(20); + DoubleArray InterfaceSpeed(20); + DoubleArray NormalVector(60); + DoubleArray van(3); + DoubleArray vaw(3); + DoubleArray vawn(3); + DoubleArray Gwn(6); + DoubleArray Gns(6); + DoubleArray Gws(6); + + int c; + //........................................................................... + int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo + IntArray cubeList(3,ncubes); + pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz); + //........................................................................... + double Cx,Cy,Cz; + double dist1,dist2; + // Extra copies of phase indicator needed to compute time derivatives on CPU + DoubleArray Phase_tminus(Nx,Ny,Nz); + DoubleArray Phase_tplus(Nx,Ny,Nz); + DoubleArray dPdt(Nx,Ny,Nz); + Cx = Cy = Cz = Nz*0.5; + for (k=0; k 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 ){ + // volume the excludes the interfacial region + vol_n += 1.0; + // pressure + pan += Press(i,j,k); + // velocity + van(0) += Vel_x(i,j,k); + van(1) += Vel_y(i,j,k); + van(2) += Vel_z(i,j,k); + } + + // volume averages over the wetting phase + if ( Phase(i,j,k) < -0.9999 ){ + // volume the excludes the interfacial region + vol_w += 1.0; + // pressure + paw += Press(i,j,k); + // velocity + vaw(0) += Vel_x(i,j,k); + vaw(1) += Vel_y(i,j,k); + vaw(2) += Vel_z(i,j,k); + } + } + } + } + } + + // End of the loop to set the values + awn = aws = ans = lwns = 0.0; + nwp_volume = 0.0; + As = 0.0; + // Compute phase averages + pan = paw = 0.0; + vaw(0) = vaw(1) = vaw(2) = 0.0; + Gwn(0) = Gwn(1) = Gwn(2) = Gwn(3) = Gwn(4) = Gwn(5) = 0.0; + Gws(0) = Gws(1) = Gws(2) = Gws(3) = Gws(4) = Gws(5) = 0.0; + Gns(0) = Gns(1) = Gns(2) = Gns(3) = Gns(4) = Gns(5) = 0.0; + + vol_w = vol_n =0.0; + + for (c=0;c