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

Interface speed and orientation tensor validated

Browse Source
This commit is contained in:
James McClure
2013-12-14 12:23:41 -05:00
parent cefcd9a5b9
commit a7ce65eb29
3 changed files with 506 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

172
gpu/lb2_Color_wia_mpi.cpp
View File

@@ -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<count; idx++){
n = list[idx];
sendbuf[idx] = Values.data[n];
}
}
inline void UnpackMeshData(int *list, int count, double *recvbuf, DoubleArray &Values){
// Fill in the phase ID values from neighboring processors
// This unpacks the values once they have been recieved from neighbors
int idx,n;
for (idx=0; idx<count; idx++){
n = list[idx];
Values.data[n] = recvbuf[idx];
}
}
//***************************************************************************************
inline void CommunicateMeshHalo()
{
sendtag = recvtag = 7;
PackID(sendList_x, sendCount_x ,sendID_x, id);
PackID(sendList_X, sendCount_X ,sendID_X, id);
PackID(sendList_y, sendCount_y ,sendID_y, id);
PackID(sendList_Y, sendCount_Y ,sendID_Y, id);
PackID(sendList_z, sendCount_z ,sendID_z, id);
PackID(sendList_Z, sendCount_Z ,sendID_Z, id);
PackID(sendList_xy, sendCount_xy ,sendID_xy, id);
PackID(sendList_Xy, sendCount_Xy ,sendID_Xy, id);
PackID(sendList_xY, sendCount_xY ,sendID_xY, id);
PackID(sendList_XY, sendCount_XY ,sendID_XY, id);
PackID(sendList_xz, sendCount_xz ,sendID_xz, id);
PackID(sendList_Xz, sendCount_Xz ,sendID_Xz, id);
PackID(sendList_xZ, sendCount_xZ ,sendID_xZ, id);
PackID(sendList_XZ, sendCount_XZ ,sendID_XZ, id);
PackID(sendList_yz, sendCount_yz ,sendID_yz, id);
PackID(sendList_Yz, sendCount_Yz ,sendID_Yz, id);
PackID(sendList_yZ, sendCount_yZ ,sendID_yZ, id);
PackID(sendList_YZ, sendCount_YZ ,sendID_YZ, id);
//......................................................................................
MPI_Sendrecv(sendID_x,sendCount_x,MPI_CHAR,rank_x,sendtag,
recvID_X,recvCount_X,MPI_CHAR,rank_X,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_X,sendCount_X,MPI_CHAR,rank_X,sendtag,
recvID_x,recvCount_x,MPI_CHAR,rank_x,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_y,sendCount_y,MPI_CHAR,rank_y,sendtag,
recvID_Y,recvCount_Y,MPI_CHAR,rank_Y,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Y,sendCount_Y,MPI_CHAR,rank_Y,sendtag,
recvID_y,recvCount_y,MPI_CHAR,rank_y,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_z,sendCount_z,MPI_CHAR,rank_z,sendtag,
recvID_Z,recvCount_Z,MPI_CHAR,rank_Z,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Z,sendCount_Z,MPI_CHAR,rank_Z,sendtag,
recvID_z,recvCount_z,MPI_CHAR,rank_z,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xy,sendCount_xy,MPI_CHAR,rank_xy,sendtag,
recvID_XY,recvCount_XY,MPI_CHAR,rank_XY,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_XY,sendCount_XY,MPI_CHAR,rank_XY,sendtag,
recvID_xy,recvCount_xy,MPI_CHAR,rank_xy,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Xy,sendCount_Xy,MPI_CHAR,rank_Xy,sendtag,
recvID_xY,recvCount_xY,MPI_CHAR,rank_xY,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xY,sendCount_xY,MPI_CHAR,rank_xY,sendtag,
recvID_Xy,recvCount_Xy,MPI_CHAR,rank_Xy,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xz,sendCount_xz,MPI_CHAR,rank_xz,sendtag,
recvID_XZ,recvCount_XZ,MPI_CHAR,rank_XZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_XZ,sendCount_XZ,MPI_CHAR,rank_XZ,sendtag,
recvID_xz,recvCount_xz,MPI_CHAR,rank_xz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Xz,sendCount_Xz,MPI_CHAR,rank_Xz,sendtag,
recvID_xZ,recvCount_xZ,MPI_CHAR,rank_xZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_xZ,sendCount_xZ,MPI_CHAR,rank_xZ,sendtag,
recvID_Xz,recvCount_Xz,MPI_CHAR,rank_Xz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_yz,sendCount_yz,MPI_CHAR,rank_yz,sendtag,
recvID_YZ,recvCount_YZ,MPI_CHAR,rank_YZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_YZ,sendCount_YZ,MPI_CHAR,rank_YZ,sendtag,
recvID_yz,recvCount_yz,MPI_CHAR,rank_yz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_Yz,sendCount_Yz,MPI_CHAR,rank_Yz,sendtag,
recvID_yZ,recvCount_yZ,MPI_CHAR,rank_yZ,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendID_yZ,sendCount_yZ,MPI_CHAR,rank_yZ,sendtag,
recvID_Yz,recvCount_Yz,MPI_CHAR,rank_Yz,recvtag,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
//......................................................................................
UnpackID(recvList_x, recvCount_x ,recvID_x, id);
UnpackID(recvList_X, recvCount_X ,recvID_X, id);
UnpackID(recvList_y, recvCount_y ,recvID_y, id);
UnpackID(recvList_Y, recvCount_Y ,recvID_Y, id);
UnpackID(recvList_z, recvCount_z ,recvID_z, id);
UnpackID(recvList_Z, recvCount_Z ,recvID_Z, id);
UnpackID(recvList_xy, recvCount_xy ,recvID_xy, id);
UnpackID(recvList_Xy, recvCount_Xy ,recvID_Xy, id);
UnpackID(recvList_xY, recvCount_xY ,recvID_xY, id);
UnpackID(recvList_XY, recvCount_XY ,recvID_XY, id);
UnpackID(recvList_xz, recvCount_xz ,recvID_xz, id);
UnpackID(recvList_Xz, recvCount_Xz ,recvID_Xz, id);
UnpackID(recvList_xZ, recvCount_xZ ,recvID_xZ, id);
UnpackID(recvList_XZ, recvCount_XZ ,recvID_XZ, id);
UnpackID(recvList_yz, recvCount_yz ,recvID_yz, id);
UnpackID(recvList_Yz, recvCount_Yz ,recvID_Yz, id);
UnpackID(recvList_yZ, recvCount_yZ ,recvID_yZ, id);
UnpackID(recvList_YZ, recvCount_YZ ,recvID_YZ, id);
}
//***************************************************************************************
int main(int argc, char **argv)
{
//*****************************************
@@ -1223,18 +1326,53 @@ int main(int argc, char **argv)
UnpackID(recvList_Yz, recvCount_Yz ,recvID_Yz, id);
UnpackID(recvList_yZ, recvCount_yZ ,recvID_yZ, id);
UnpackID(recvList_YZ, recvCount_YZ ,recvID_YZ, id);
//.....................................................................................
/* // Once the ID is saved, free memory allocated to the buffers (no longer needed)
//......................................................................................
free(sendID_x); free(sendID_X); free(sendID_y); free(sendID_Y); free(sendID_z); free(sendID_Z);
free(sendID_xy); free(sendID_XY); free(sendID_xY); free(sendID_Xy);
free(sendID_xz); free(sendID_XZ); free(sendID_xZ); free(sendID_Xz);
free(sendID_yz); free(sendID_YZ); free(sendID_yZ); free(sendID_Yz);
free(recvID_x); free(recvID_X); free(recvID_y); free(recvID_Y); free(recvID_z); free(recvID_Z);
free(recvID_xy); free(recvID_XY); free(recvID_xY); free(recvID_Xy);
free(recvID_xz); free(recvID_XZ); free(recvID_xZ); free(recvID_Xz);
free(recvID_yz); free(recvID_YZ); free(recvID_yZ); free(recvID_Yz);
*/ //......................................................................................
// Fill in the phase MeshData from neighboring processors
double *sendMeshData_x, *sendMeshData_y, *sendMeshData_z, *sendMeshData_X, *sendMeshData_Y, *sendMeshData_Z;
double *sendMeshData_xy, *sendMeshData_yz, *sendMeshData_xz, *sendMeshData_Xy, *sendMeshData_Yz, *sendMeshData_xZ;
double *sendMeshData_xY, *sendMeshData_yZ, *sendMeshData_Xz, *sendMeshData_XY, *sendMeshData_YZ, *sendMeshData_XZ;
double *recvMeshData_x, *recvMeshData_y, *recvMeshData_z, *recvMeshData_X, *recvMeshData_Y, *recvMeshData_Z;
double *recvMeshData_xy, *recvMeshData_yz, *recvMeshData_xz, *recvMeshData_Xy, *recvMeshData_Yz, *recvMeshData_xZ;
double *recvMeshData_xY, *recvMeshData_yZ, *recvMeshData_Xz, *recvMeshData_XY, *recvMeshData_YZ, *recvMeshData_XZ;
// send buffers
sendMeshData_x = new double [sendCount_x];
sendMeshData_y = new double [sendCount_y];
sendMeshData_z = new double [sendCount_z];
sendMeshData_X = new double [sendCount_X];
sendMeshData_Y = new double [sendCount_Y];
sendMeshData_Z = new double [sendCount_Z];
sendMeshData_xy = new double [sendCount_xy];
sendMeshData_yz = new double [sendCount_yz];
sendMeshData_xz = new double [sendCount_xz];
sendMeshData_Xy = new double [sendCount_Xy];
sendMeshData_Yz = new double [sendCount_Yz];
sendMeshData_xZ = new double [sendCount_xZ];
sendMeshData_xY = new double [sendCount_xY];
sendMeshData_yZ = new double [sendCount_yZ];
sendMeshData_Xz = new double [sendCount_Xz];
sendMeshData_XY = new double [sendCount_XY];
sendMeshData_YZ = new double [sendCount_YZ];
sendMeshData_XZ = new double [sendCount_XZ];
//......................................................................................
// recv buffers
recvMeshData_x = new double [recvCount_x];
recvMeshData_y = new double [recvCount_y];
recvMeshData_z = new double [recvCount_z];
recvMeshData_X = new double [recvCount_X];
recvMeshData_Y = new double [recvCount_Y];
recvMeshData_Z = new double [recvCount_Z];
recvMeshData_xy = new double [recvCount_xy];
recvMeshData_yz = new double [recvCount_yz];
recvMeshData_xz = new double [recvCount_xz];
recvMeshData_Xy = new double [recvCount_Xy];
recvMeshData_xZ = new double [recvCount_xZ];
recvMeshData_xY = new double [recvCount_xY];
recvMeshData_yZ = new double [recvCount_yZ];
recvMeshData_Yz = new double [recvCount_Yz];
recvMeshData_Xz = new double [recvCount_Xz];
recvMeshData_XY = new double [recvCount_XY];
recvMeshData_YZ = new double [recvCount_YZ];
recvMeshData_XZ = new double [recvCount_XZ];
if (rank==0) printf ("Devices are ready to communicate. \n");
MPI_Barrier(MPI_COMM_WORLD);
@@ -1342,6 +1480,9 @@ int main(int argc, char **argv)
IntArray nws_seg(2,20);
DTMutableList<Point> 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<N; n++) dPdt(n) = 0.1*(Phase_plus(n) - Phase_tminus(n));
//...........................................................................
// 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);
@@ -1929,6 +2071,8 @@ int main(int argc, char **argv)
// Compute the mesh curvature of the phase indicator field
pmmc_MeshCurvature(Phase, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
//...........................................................................
// Fill in the halo region for the mesh gradients and curvature
//...........................................................................
// Compute areas using porous medium marching cubes algorithm
@@ -1938,7 +2082,6 @@ int main(int argc, char **argv)
nwp_volume = 0.0;
As = 0.0;
// Compute phase averages
p_n = p_w = 0.0;
vx_w = vy_w = vz_w = 0.0;
@@ -2005,6 +2148,9 @@ int main(int argc, char **argv)
// Integrate the mean curvature
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
//...........................................................................
// Compute the Interfacial Areas, Common Line length
awn += pmmc_CubeSurfaceArea(nw_pts,nw_tris,n_nw_tris);

16
include/pmmc.h
View File

@@ -3792,7 +3792,7 @@ inline double pmmc_CubeSurfaceOrientation(DoubleArray &Orientation, DTMutableLis
{
int r;
double temp,area,s,s1,s2,s3;
double nx,ny,nz,norm;
double nx,ny,nz,normsq;
Point A,B,C;
area = 0.0;
for (r=0;r<ntris;r++){
@@ -3803,7 +3803,7 @@ inline double pmmc_CubeSurfaceOrientation(DoubleArray &Orientation, DTMutableLis
nx = (B.y-A.y)*(C.z-A.z) - (B.z-A.z)*(C.y-A.y);
ny = (B.z-A.z)*(C.x-A.x) - (B.x-A.x)*(C.z-A.z);
nz = (B.x-A.x)*(C.y-A.y) - (B.y-A.y)*(C.x-A.x);
norm = sqrt(nx*nx+ny*ny+nz*nz);
normsq = 1.0/(nx*nx+ny*ny+nz*nz);
// Compute length of sides (assume dx=dy=dz)
s1 = 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));
s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
@@ -3813,12 +3813,12 @@ inline double pmmc_CubeSurfaceOrientation(DoubleArray &Orientation, DTMutableLis
if (temp > 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;

337
pmmc/Analysis.cpp Normal file
View File

@@ -0,0 +1,337 @@
#include <iostream>
#include <math.h>
#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<Point> nw_pts(20);
DTMutableList<Point> ns_pts(20);
DTMutableList<Point> ws_pts(20);
DTMutableList<Point> 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<Point> tmp(20);
// Initialize arrays for local solid surface
DTMutableList<Point> 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<Point> 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<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD;
Phase_tminus(i,j,k) = dist2;
}
}
}
Cz += SPEED;
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
dist1 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)) - RADIUS;
dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD;
SignDist(i,j,k) = -dist1;
Phase(i,j,k) = dist2;
}
}
}
Cz += SPEED;
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD;
Phase_tplus(i,j,k) = dist2;
}
}
}
//...........................................................................
// 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);
pmmc_MeshGradient(SignDist,Sx,Sy,Sz,Nx,Ny,Nz);
pmmc_MeshCurvature(Phase, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
/// Compute volume averages
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
if ( SignDist(i,j,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<ncubes;c++){
// Get cube from the list
i = cubeList(0,c);
j = cubeList(1,c);
k = cubeList(2,c);
// Run PMMC
n_local_sol_tris = 0;
n_local_sol_pts = 0;
n_local_nws_pts = 0;
n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(SignDist, Phase, solid_isovalue, fluid_isovalue,
nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
i, j, k, Nx, Ny, Nz);
// Compute the velocity of the wn interface
pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
// 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);
// 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);
// Compute the surface orientation and the interfacial area
awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
//*******************************************************************
// Compute the Interfacial Areas, Common Line length
// awn += pmmc_CubeSurfaceArea(nw_pts,nw_tris,n_nw_tris);
// ans += pmmc_CubeSurfaceArea(ns_pts,ns_tris,n_ns_tris);
// aws += pmmc_CubeSurfaceArea(ws_pts,ws_tris,n_ws_tris);
As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
}
Jwn /= awn;
efawns /= lwns;
for (i=0; i<3; i++) vawn(i) /= awn;
for (i=0; i<6; i++) Gwn(i) /= awn;
for (i=0; i<6; i++) Gns(i) /= ans;
for (i=0; i<6; i++) Gws(i) /= aws;
printf("--------------------------------------------------------------------------------------\n");
printf("sw pw pn vw[x, y, z] vn[x, y, z] "); // Volume averages
printf("awn ans aws Jwn vwn[x, y, z] lwns efawns "); // Interface and common curve averages
printf("Gwn [xx, yy, zz, xy, xz, yz] "); // Orientation tensors
printf("Gws [xx, yy, zz, xy, xz, yz] ");
printf("Gns [xx, yy, zz, xy, xz, yz] \n");
printf("--------------------------------------------------------------------------------------\n");
printf("%.5g %.5g %.5g ",sw,paw,pan); // saturation and pressure
printf("%.5g %.5g %.5g ",vaw(0),vaw(1),vaw(2)); // average velocity of w phase
printf("%.5g %.5g %.5g ",van(0),van(1),van(2)); // average velocity of n phase
printf("%.5g %.5g %.5g ",awn,ans,aws); // interfacial areas
printf("%.5g ",Jwn); // curvature of wn interface
printf("%.5g %.5g %.5g ",vawn(0),vawn(1),vawn(2)); // velocity of wn interface
printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
Gwn(0),Gwn(1),Gwn(2),Gwn(3),Gwn(4),Gwn(5)); // orientation of wn interface
printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
Gns(0),Gns(1),Gns(2),Gns(3),Gns(4),Gns(5)); // orientation of ns interface
printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
Gws(0),Gws(1),Gws(2),Gws(3),Gws(4),Gws(5)); // orientation of ws interface
/* printf("-------------------------------- \n");
printf("NWP volume = %f \n", nwp_volume);
printf("Area wn = %f, Analytical = %f \n", awn,2*PI*RADIUS*RADIUS);
printf("Area ns = %f, Analytical = %f \n", ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT);
printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT);
printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2));
printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS);
// printf("Cos(theta_wns) = %f, Analytical = %f \n",efawns/lwns,1.0*RADIUS/CAPRAD);
printf("Interface Velocity = %f,%f,%f \n",vawn(0)/awn,vawn(1)/awn,vawn(2)/awn);
printf("-------------------------------- \n");
//.........................................................................
FILE *PHASE;
PHASE = fopen("Phase.in","wb");
fwrite(Phase,8,SIZE,PHASE);
fclose(PHASE);
FILE *SOLID;
SOLID = fopen("Distance.in","wb");
fwrite(Solid,8,SIZE,SOLID);
fclose(SOLID);
*/
}