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

Merge branch 'master' into cmake

Browse Source
This commit is contained in:
Mark Berrill 2014-01-28 13:29:55 -05:00
commit df295faa84
21 changed files with 7698 additions and 525 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

833
cpu/Color.cpp
View File

File diff suppressed because it is too large Load Diff

41
cpu/Color.h
View File

@ -1,21 +1,20 @@
extern void InitDenColor(char *ID, double *Den, double *Phi, double das, double dbs, int N);
extern void InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz);
extern void Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz);
//*************************************************************************
//*************************************************************************
extern void PressureBC_inlet(double *disteven, double *distodd, double din,
int Nx, int Ny, int Nz);
extern void PressureBC_outlet(double *disteven, double *distodd, double dout,
int Nx, int Ny, int Nz, int S, int outlet);
//*************************************************************************
extern void ComputeColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz);
//*************************************************************************
extern void ColorCollide( char *ID, double *disteven, double *distodd, double *ColorGrad,
double *Velocity, int Nx, int Ny, int Nz, double rlx_setA, double rlx_setB,
double alpha, double beta, double Fx, double Fy, double Fz, bool pBC);
//*************************************************************************
extern void DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity,
double beta, int Nx, int Ny, int Nz, bool pBC);
extern void ComputePhi(char *ID, double *Phi, double *Copy, double *Den, int N);
extern "C" void dvc_InitDenColor(char *ID, double *Den, double *Phi, double das, double dbs, int Nx, int Ny, int Nz, int S);
extern "C" void dvc_InitDenColorDistance(char *ID, double *Den, double *Phi, double *Distance,
double das, double dbs, double beta, double xp, int Nx, int Ny, int Nz, int S);
extern "C" void dvc_Compute_VELOCITY(char *ID, double *disteven, double *distodd, double *vel, int Nx, int Ny, int Nz, int S);
extern "C" void dvc_ComputePressureD3Q19(char *ID, double *disteven, double *distodd, double *Pressure,
int Nx, int Ny, int Nz, int S);
extern "C" void dvc_PressureBC_inlet(double *disteven, double *distodd, double din,
int Nx, int Ny, int Nz, int S);
extern "C" void dvc_PressureBC_outlet(double *disteven, double *distodd, double dout,
int Nx, int Ny, int Nz, int S, int outlet);
extern "C" void dvc_ComputeColorGradient(char *ID, double *phi, double *ColorGrad, int Nx, int Ny, int Nz, int S);
extern "C" void dvc_ColorCollide( char *ID, double *disteven, double *distodd, double *ColorGrad,
double *Velocity, int Nx, int Ny, int Nz, int S,double rlx_setA, double rlx_setB,
double alpha, double beta, double Fx, double Fy, double Fz, bool pBC);
extern "C" void dvc_ColorCollideOpt( char *ID, double *disteven, double *distodd, double *phi, double *ColorGrad,
double *Velocity, int Nx, int Ny, int Nz, int S,double rlx_setA, double rlx_setB,
double alpha, double beta, double Fx, double Fy, double Fz);
extern "C" void dvc_DensityStreamD3Q7(char *ID, double *Den, double *Copy, double *Phi, double *ColorGrad, double *Velocity,
double beta, int Nx, int Ny, int Nz, bool pBC, int S);
extern "C" void dvc_ComputePhi(char *ID, double *Phi, double *Copy, double *Den, int N, int S);

46
cpu/D3Q19.cpp
View File

@ -1,4 +1,4 @@
extern void PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){
extern "C" void dvc_PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N){
//....................................................................................
// Pack distribution q into the send buffer for the listed lattice sites
// dist may be even or odd distributions stored by stream layout
@ -10,9 +10,7 @@ extern void PackDist(int q, int *list, int start, int count, double *sendbuf, do
}
}
extern void MapRecvDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count,
extern "C" void dvc_UnpackDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count,
double *recvbuf, double *dist, int Nx, int Ny, int Nz){
//....................................................................................
// Unack distribution from the recv buffer
@ -54,8 +52,46 @@ extern void MapRecvDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start,
}
}
extern "C" void dvc_InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz, int S)
{
int n,N;
N = Nx*Ny*Nz;
for (n=0; n<N; n++){
if (ID[n] > 0){
f_even[n] = 0.3333333333333333;
f_odd[n] = 0.055555555555555555; //double(100*n)+1.f;
f_even[N+n] = 0.055555555555555555; //double(100*n)+2.f;
f_odd[N+n] = 0.055555555555555555; //double(100*n)+3.f;
f_even[2*N+n] = 0.055555555555555555; //double(100*n)+4.f;
f_odd[2*N+n] = 0.055555555555555555; //double(100*n)+5.f;
f_even[3*N+n] = 0.055555555555555555; //double(100*n)+6.f;
f_odd[3*N+n] = 0.0277777777777778; //double(100*n)+7.f;
f_even[4*N+n] = 0.0277777777777778; //double(100*n)+8.f;
f_odd[4*N+n] = 0.0277777777777778; //double(100*n)+9.f;
f_even[5*N+n] = 0.0277777777777778; //double(100*n)+10.f;
f_odd[5*N+n] = 0.0277777777777778; //double(100*n)+11.f;
f_even[6*N+n] = 0.0277777777777778; //double(100*n)+12.f;
f_odd[6*N+n] = 0.0277777777777778; //double(100*n)+13.f;
f_even[7*N+n] = 0.0277777777777778; //double(100*n)+14.f;
f_odd[7*N+n] = 0.0277777777777778; //double(100*n)+15.f;
f_even[8*N+n] = 0.0277777777777778; //double(100*n)+16.f;
f_odd[8*N+n] = 0.0277777777777778; //double(100*n)+17.f;
f_even[9*N+n] = 0.0277777777777778; //double(100*n)+18.f;
}
else{
for(int q=0; q<9; q++){
f_even[q*N+n] = -1.0;
f_odd[q*N+n] = -1.0;
}
f_even[9*N+n] = -1.0;
}
}
}
//*************************************************************************
extern void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz)
extern "C" void dvc_SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz, int S)
{
int n,nn,N;
// distributions

8
cpu/D3Q19.h
View File

@ -1,6 +1,6 @@
extern void PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N);
extern void MapRecvDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count,
extern "C" void dvc_PackDist(int q, int *list, int start, int count, double *sendbuf, double *dist, int N);
extern "C" void dvc_UnpackDist(int q, int Cqx, int Cqy, int Cqz, int *list, int start, int count,
double *recvbuf, double *dist, int Nx, int Ny, int Nz);
//*************************************************************************
extern void SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz);
extern "C" void dvc_InitD3Q19(char *ID, double *f_even, double *f_odd, int Nx, int Ny, int Nz, int S);
extern "C" void dvc_SwapD3Q19(char *ID, double *disteven, double *distodd, int Nx, int Ny, int Nz, int S);

10
cpu/D3Q7.cpp
View File

@ -1,6 +1,6 @@
// GPU Functions for D3Q7 Lattice Boltzmann Methods
// CPU Functions for D3Q7 Lattice Boltzmann Methods
extern void PackValues(int *list, int count, double *sendbuf, double *Data, int N){
extern "C" void dvc_PackValues(int *list, int count, double *sendbuf, double *Data, int N){
//....................................................................................
// Pack distribution q into the send buffer for the listed lattice sites
// dist may be even or odd distributions stored by stream layout
@ -11,7 +11,7 @@ extern void PackValues(int *list, int count, double *sendbuf, double *Data, int
sendbuf[idx] = Data[n];
}
}
extern void UnpackValues(int *list, int count, double *recvbuf, double *Data, int N){
extern "C" void dvc_UnpackValues(int *list, int count, double *recvbuf, double *Data, int N){
//....................................................................................
// Pack distribution q into the send buffer for the listed lattice sites
// dist may be even or odd distributions stored by stream layout
@ -23,7 +23,7 @@ extern void UnpackValues(int *list, int count, double *recvbuf, double *Data, in
}
}
extern void PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N){
extern "C" void dvc_PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N){
//....................................................................................
// Pack distribution into the send buffer for the listed lattice sites
//....................................................................................
@ -38,7 +38,7 @@ extern void PackDenD3Q7(int *list, int count, double *sendbuf, int number, doubl
}
extern void UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N){
extern "C" void dvc_UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N){
//....................................................................................
// Unack distribution from the recv buffer
// Sum to the existing density value

8
cpu/D3Q7.h
View File

@ -1,9 +1,9 @@
// CPU Functions for D3Q7 Lattice Boltzmann Methods
extern void PackValues(int *list, int count, double *sendbuf, double *Data, int N);
extern "C" void dvc_PackValues(int *list, int count, double *sendbuf, double *Data, int N);
extern void UnpackValues(int *list, int count, double *recvbuf, double *Data, int N);
extern "C" void dvc_UnpackValues(int *list, int count, double *recvbuf, double *Data, int N);
extern void PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N);
extern "C" void dvc_PackDenD3Q7(int *list, int count, double *sendbuf, int number, double *Data, int N);
extern void UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N);
extern "C" void dvc_UnpackDenD3Q7(int *list, int count, double *recvbuf, int number, double *Data, int N);

28
cpu/Extras.cpp Normal file
View File

@ -0,0 +1,28 @@
// Basic cuda functions callable from C/C++ code
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size){
//cudaMalloc(address,size);
(*address) = malloc(size);
if (*address==NULL){
printf("Memory allocation failed! \n");
}
}
extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size){
// cudaMemcpy(dest,source,size,cudaMemcpyHostToDevice);
memcpy(dest, source, size);
}
extern "C" void dvc_CopyToHost(void* dest, void* source, size_t size){
// cudaMemcpy(dest,source,size,cudaMemcpyDeviceToHost);
memcpy(dest, source, size);
}
extern "C" void dvc_Barrier(){
// cudaDeviceSynchronize();
}

7
cpu/Extras.h Normal file
View File

@ -0,0 +1,7 @@
extern "C" void dvc_AllocateDeviceMemory(void** address, size_t size);
extern "C" void dvc_CopyToDevice(void* dest, void* source, size_t size);
extern "C" void dvc_CopyToHost(void* dest, void* source, size_t size);
extern "C" void dvc_Barrier();

26
cpu/Makefile
View File

@ -1,8 +1,17 @@
CXX=mpicxx
FLAGS=-O3
INC=../include
ColorLBM-cpu:D3Q19.o D3Q7.o Color.o lb2_Color_mpi.o
$(CXX) $(FLAGS) -o ColorLBM-cpu lb2_Color_mpi.o D3Q19.o D3Q7.o Color.o
all:ColorLBM-cpu ColorLBM-CBUB
ColorLBM-cpu:Extras.o D3Q19.o D3Q7.o Color.o lb2_Color_wia_mpi.o
$(CXX) $(FLAGS) -I$(INC) -o ColorLBM-cpu lb2_Color_wia_mpi.o D3Q19.o D3Q7.o Color.o Extras.o
ColorLBM-CBUB:Extras.o D3Q19.o D3Q7.o Color.o lb2_Color_wia_mpi-CBUB.o
$(CXX) $(FLAGS) -I$(INC) -DCBUB -o ColorLBM-CBUB lb2_Color_wia_mpi-CBUB.o D3Q19.o D3Q7.o Color.o Extras.o
ColorLBM-BUB:Extras.o D3Q19.o D3Q7.o Color.o lb2_Color_wia_mpi_bubble.o
$(CXX) $(FLAGS) -I$(INC) -DCBUB -o ColorLBM-BUB lb2_Color_wia_mpi_bubble.o D3Q19.o D3Q7.o Color.o Extras.o
D3Q19.o:D3Q19.cpp
$(CXX) $(FLAGS) -c -o D3Q19.o D3Q19.cpp
@ -10,11 +19,20 @@ D3Q19.o:D3Q19.cpp
D3Q7.o:D3Q7.cpp
$(CXX) $(FLAGS) -c -o D3Q7.o D3Q7.cpp
Extras.o:Extras.cpp
$(CXX) $(FLAGS) -c -o Extras.o Extras.cpp
Color.o:Color.cpp
$(CXX) $(FLAGS) -c -o Color.o Color.cpp
lb2_Color_mpi.o:lb2_Color_mpi.cpp
$(CXX) $(FLAGS) -c -o lb2_Color_mpi.o lb2_Color_mpi.cpp
lb2_Color_wia_mpi.o:lb2_Color_wia_mpi.cpp
$(CXX) $(FLAGS) -I$(INC) -c -o lb2_Color_wia_mpi.o lb2_Color_wia_mpi.cpp
lb2_Color_wia_mpi-CBUB.o:lb2_Color_wia_mpi.cpp
$(CXX) $(FLAGS) -I$(INC) -DCBUB -c -o lb2_Color_wia_mpi-CBUB.o lb2_Color_wia_mpi.cpp
lb2_Color_wia_mpi_bubble.o:lb2_Color_wia_mpi_bubble.cpp
$(CXX) $(FLAGS) -I$(INC) -DCBUB -c -o lb2_Color_wia_mpi_bubble.o lb2_Color_wia_mpi_bubble.cpp
#MRT-MPI.o:lb1_MRT_mpi.cpp
# $(CXX) -c -o MRT-MPI.o lb1_MRT_mpi.cpp

2447
cpu/lb2_Color_wia_mpi.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

2492
cpu/lb2_Color_wia_mpi_bubble.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

783
gpu/exe/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,136 @@ 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(DoubleArray &MeshData, MPI_Comm Communicator,
double *sendbuf_x,double *sendbuf_y,double *sendbuf_z,double *sendbuf_X,double *sendbuf_Y,double *sendbuf_Z,
double *sendbuf_xy,double *sendbuf_XY,double *sendbuf_xY,double *sendbuf_Xy,
double *sendbuf_xz,double *sendbuf_XZ,double *sendbuf_xZ,double *sendbuf_Xz,
double *sendbuf_yz,double *sendbuf_YZ,double *sendbuf_yZ,double *sendbuf_Yz,
double *recvbuf_x,double *recvbuf_y,double *recvbuf_z,double *recvbuf_X,double *recvbuf_Y,double *recvbuf_Z,
double *recvbuf_xy,double *recvbuf_XY,double *recvbuf_xY,double *recvbuf_Xy,
double *recvbuf_xz,double *recvbuf_XZ,double *recvbuf_xZ,double *recvbuf_Xz,
double *recvbuf_yz,double *recvbuf_YZ,double *recvbuf_yZ,double *recvbuf_Yz,
int *sendList_x,int *sendList_y,int *sendList_z,int *sendList_X,int *sendList_Y,int *sendList_Z,
int *sendList_xy,int *sendList_XY,int *sendList_xY,int *sendList_Xy,
int *sendList_xz,int *sendList_XZ,int *sendList_xZ,int *sendList_Xz,
int *sendList_yz,int *sendList_YZ,int *sendList_yZ,int *sendList_Yz,
int sendCount_x,int sendCount_y,int sendCount_z,int sendCount_X,int sendCount_Y,int sendCount_Z,
int sendCount_xy,int sendCount_XY,int sendCount_xY,int sendCount_Xy,
int sendCount_xz,int sendCount_XZ,int sendCount_xZ,int sendCount_Xz,
int sendCount_yz,int sendCount_YZ,int sendCount_yZ,int sendCount_Yz,
int *recvList_x,int *recvList_y,int *recvList_z,int *recvList_X,int *recvList_Y,int *recvList_Z,
int *recvList_xy,int *recvList_XY,int *recvList_xY,int *recvList_Xy,
int *recvList_xz,int *recvList_XZ,int *recvList_xZ,int *recvList_Xz,
int *recvList_yz,int *recvList_YZ,int *recvList_yZ,int *recvList_Yz,
int recvCount_x,int recvCount_y,int recvCount_z,int recvCount_X,int recvCount_Y,int recvCount_Z,
int recvCount_xy,int recvCount_XY,int recvCount_xY,int recvCount_Xy,
int recvCount_xz,int recvCount_XZ,int recvCount_xZ,int recvCount_Xz,
int recvCount_yz,int recvCount_YZ,int recvCount_yZ,int recvCount_Yz,
int rank_x,int rank_y,int rank_z,int rank_X,int rank_Y,int rank_Z,int rank_xy,int rank_XY,int rank_xY,
int rank_Xy,int rank_xz,int rank_XZ,int rank_xZ,int rank_Xz,int rank_yz,int rank_YZ,int rank_yZ,int rank_Yz)
{
int sendtag, recvtag;
sendtag = recvtag = 7;
PackMeshData(sendList_x, sendCount_x ,sendbuf_x, MeshData);
PackMeshData(sendList_X, sendCount_X ,sendbuf_X, MeshData);
PackMeshData(sendList_y, sendCount_y ,sendbuf_y, MeshData);
PackMeshData(sendList_Y, sendCount_Y ,sendbuf_Y, MeshData);
PackMeshData(sendList_z, sendCount_z ,sendbuf_z, MeshData);
PackMeshData(sendList_Z, sendCount_Z ,sendbuf_Z, MeshData);
PackMeshData(sendList_xy, sendCount_xy ,sendbuf_xy, MeshData);
PackMeshData(sendList_Xy, sendCount_Xy ,sendbuf_Xy, MeshData);
PackMeshData(sendList_xY, sendCount_xY ,sendbuf_xY, MeshData);
PackMeshData(sendList_XY, sendCount_XY ,sendbuf_XY, MeshData);
PackMeshData(sendList_xz, sendCount_xz ,sendbuf_xz, MeshData);
PackMeshData(sendList_Xz, sendCount_Xz ,sendbuf_Xz, MeshData);
PackMeshData(sendList_xZ, sendCount_xZ ,sendbuf_xZ, MeshData);
PackMeshData(sendList_XZ, sendCount_XZ ,sendbuf_XZ, MeshData);
PackMeshData(sendList_yz, sendCount_yz ,sendbuf_yz, MeshData);
PackMeshData(sendList_Yz, sendCount_Yz ,sendbuf_Yz, MeshData);
PackMeshData(sendList_yZ, sendCount_yZ ,sendbuf_yZ, MeshData);
PackMeshData(sendList_YZ, sendCount_YZ ,sendbuf_YZ, MeshData);
//......................................................................................
MPI_Sendrecv(sendbuf_x,sendCount_x,MPI_CHAR,rank_x,sendtag,
recvbuf_X,recvCount_X,MPI_CHAR,rank_X,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_X,sendCount_X,MPI_CHAR,rank_X,sendtag,
recvbuf_x,recvCount_x,MPI_CHAR,rank_x,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_y,sendCount_y,MPI_CHAR,rank_y,sendtag,
recvbuf_Y,recvCount_Y,MPI_CHAR,rank_Y,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Y,sendCount_Y,MPI_CHAR,rank_Y,sendtag,
recvbuf_y,recvCount_y,MPI_CHAR,rank_y,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_z,sendCount_z,MPI_CHAR,rank_z,sendtag,
recvbuf_Z,recvCount_Z,MPI_CHAR,rank_Z,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Z,sendCount_Z,MPI_CHAR,rank_Z,sendtag,
recvbuf_z,recvCount_z,MPI_CHAR,rank_z,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xy,sendCount_xy,MPI_CHAR,rank_xy,sendtag,
recvbuf_XY,recvCount_XY,MPI_CHAR,rank_XY,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_XY,sendCount_XY,MPI_CHAR,rank_XY,sendtag,
recvbuf_xy,recvCount_xy,MPI_CHAR,rank_xy,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Xy,sendCount_Xy,MPI_CHAR,rank_Xy,sendtag,
recvbuf_xY,recvCount_xY,MPI_CHAR,rank_xY,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xY,sendCount_xY,MPI_CHAR,rank_xY,sendtag,
recvbuf_Xy,recvCount_Xy,MPI_CHAR,rank_Xy,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xz,sendCount_xz,MPI_CHAR,rank_xz,sendtag,
recvbuf_XZ,recvCount_XZ,MPI_CHAR,rank_XZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_XZ,sendCount_XZ,MPI_CHAR,rank_XZ,sendtag,
recvbuf_xz,recvCount_xz,MPI_CHAR,rank_xz,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Xz,sendCount_Xz,MPI_CHAR,rank_Xz,sendtag,
recvbuf_xZ,recvCount_xZ,MPI_CHAR,rank_xZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_xZ,sendCount_xZ,MPI_CHAR,rank_xZ,sendtag,
recvbuf_Xz,recvCount_Xz,MPI_CHAR,rank_Xz,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_yz,sendCount_yz,MPI_CHAR,rank_yz,sendtag,
recvbuf_YZ,recvCount_YZ,MPI_CHAR,rank_YZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_YZ,sendCount_YZ,MPI_CHAR,rank_YZ,sendtag,
recvbuf_yz,recvCount_yz,MPI_CHAR,rank_yz,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_Yz,sendCount_Yz,MPI_CHAR,rank_Yz,sendtag,
recvbuf_yZ,recvCount_yZ,MPI_CHAR,rank_yZ,recvtag,Communicator,MPI_STATUS_IGNORE);
MPI_Sendrecv(sendbuf_yZ,sendCount_yZ,MPI_CHAR,rank_yZ,sendtag,
recvbuf_Yz,recvCount_Yz,MPI_CHAR,rank_Yz,recvtag,Communicator,MPI_STATUS_IGNORE);
//........................................................................................
UnpackMeshData(recvList_x, recvCount_x ,recvbuf_x, MeshData);
UnpackMeshData(recvList_X, recvCount_X ,recvbuf_X, MeshData);
UnpackMeshData(recvList_y, recvCount_y ,recvbuf_y, MeshData);
UnpackMeshData(recvList_Y, recvCount_Y ,recvbuf_Y, MeshData);
UnpackMeshData(recvList_z, recvCount_z ,recvbuf_z, MeshData);
UnpackMeshData(recvList_Z, recvCount_Z ,recvbuf_Z, MeshData);
UnpackMeshData(recvList_xy, recvCount_xy ,recvbuf_xy, MeshData);
UnpackMeshData(recvList_Xy, recvCount_Xy ,recvbuf_Xy, MeshData);
UnpackMeshData(recvList_xY, recvCount_xY ,recvbuf_xY, MeshData);
UnpackMeshData(recvList_XY, recvCount_XY ,recvbuf_XY, MeshData);
UnpackMeshData(recvList_xz, recvCount_xz ,recvbuf_xz, MeshData);
UnpackMeshData(recvList_Xz, recvCount_Xz ,recvbuf_Xz, MeshData);
UnpackMeshData(recvList_xZ, recvCount_xZ ,recvbuf_xZ, MeshData);
UnpackMeshData(recvList_XZ, recvCount_XZ ,recvbuf_XZ, MeshData);
UnpackMeshData(recvList_yz, recvCount_yz ,recvbuf_yz, MeshData);
UnpackMeshData(recvList_Yz, recvCount_Yz ,recvbuf_Yz, MeshData);
UnpackMeshData(recvList_yZ, recvCount_yZ ,recvbuf_yZ, MeshData);
UnpackMeshData(recvList_YZ, recvCount_YZ ,recvbuf_YZ, MeshData);
}
//***************************************************************************************
int main(int argc, char **argv)
{
//*****************************************
@ -152,7 +281,7 @@ int main(int argc, char **argv)
double din,dout;
double wp_saturation;
bool pBC,Restart;
int i,j,k,p,q,r,n;
int i,j,k,p,n;
// pmmc threshold values
double fluid_isovalue,solid_isovalue;
@ -574,7 +703,7 @@ int main(int argc, char **argv)
#ifdef CBUB
// Initializes a constrained bubble test
double BubbleBot = 20.0; // How big to make the NWP bubble
double BubbleTop =60.0; // How big to make the NWP bubble
double BubbleTop = 60.0; // How big to make the NWP bubble
double TubeRadius = 15.0; // Radius of the capillary tube
sum=0;
for (k=0;k<Nz;k++){
@ -1223,18 +1352,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);
@ -1271,9 +1435,9 @@ int main(int argc, char **argv)
DoubleArray Phase(Nx,Ny,Nz);
// Extra copies of phi needed to compute time derivatives on CPU
double *Phi_tminus,*Phi_tplus;
Phi_tminus = new double[N];
Phi_tplus = new double[N];
DoubleArray Phase_tminus(Nx,Ny,Nz);
DoubleArray Phase_tplus(Nx,Ny,Nz);
DoubleArray dPdt(Nx,Ny,Nz);
//copies of data needed to perform checkpointing from cpu
double *cDen, *cDistEven, *cDistOdd;
@ -1282,14 +1446,19 @@ int main(int argc, char **argv)
cDistOdd = new double[9*N];
// data needed to perform CPU-based averaging
double *Vel,*Press,*Norm_x,*Norm_y,*Norm_z,*Curvature,*dphidt;
double *Vel,*Press;
Vel = new double[3*N]; // fluid velocity
Press = new double[N]; // fluid pressure
dphidt = new double[N]; // d phi / dt
Norm_x = new double[N]; // normal in the x direction
Norm_y = new double[N]; // normal in the y direction
Norm_z = new double[N]; // normal in the z direction
Curvature = new double[N]; // curvature of phi
DoubleArray MeanCurvature(Nx,Ny,Nz);
DoubleArray GaussCurvature(Nx,Ny,Nz);
DoubleArray SignDist_x(Nx,Ny,Nz); // Gradient of the signed distance
DoubleArray SignDist_y(Nx,Ny,Nz);
DoubleArray SignDist_z(Nx,Ny,Nz);
DoubleArray Phase_x(Nx,Ny,Nz); // Gradient of the phase indicator field
DoubleArray Phase_y(Nx,Ny,Nz);
DoubleArray Phase_z(Nx,Ny,Nz);
/*****************************************************************
VARIABLES FOR THE PMMC ALGORITHM
****************************************************************** */
@ -1298,25 +1467,42 @@ int main(int argc, char **argv)
//...........................................................................
// local averages (to each MPI process)
double awn,ans,aws,lwns,nwp_volume;
double vol_w, vol_n; // volumes the exclude the interfacial region
double As;
double vol_w, vol_n; // volumes the exclude the interfacial region
double sat_w;
double p_n,p_w; // local phase averaged pressure
double vx_w,vy_w,vz_w,vx_n,vy_n,vz_n; // local phase averaged velocity
double pan,paw; // local phase averaged pressure
// double vx_w,vy_w,vz_w,vx_n,vy_n,vz_n; // local phase averaged velocity
// Global averages (all processes)
double vol_w_global, vol_n_global; // volumes the exclude the interfacial region
double awn_global,ans_global,aws_global;
double lwns_global;
double efawns,efawns_global; // averaged contact angle
double Jwn,Jwn_global; // average mean curavture - wn interface
DoubleArray van(3);
DoubleArray vaw(3);
DoubleArray vawn(3);
DoubleArray Gwn(6);
DoubleArray Gns(6);
DoubleArray Gws(6);
double nwp_volume_global; // volume for the wetting phase (for saturation)
double p_n_global,p_w_global; // global phase averaged pressure
double vx_w_global,vy_w_global,vz_w_global; // global phase averaged velocity
double vx_n_global,vy_n_global,vz_n_global; // global phase averaged velocity
// double p_n_global,p_w_global; // global phase averaged pressure
// double vx_w_global,vy_w_global,vz_w_global; // global phase averaged velocity
// double vx_n_global,vy_n_global,vz_n_global; // global phase averaged velocity
double As_global;
double dEs,dAwn,dAns; // Global surface energy (calculated by rank=0)
double pan_global,paw_global; // local phase averaged pressure
DoubleArray van_global(3);
DoubleArray vaw_global(3);
DoubleArray vawn_global(3);
DoubleArray Gwn_global(6);
DoubleArray Gns_global(6);
DoubleArray Gws_global(6);
//...........................................................................
// bool add=1; // Set to false if any corners contain nw-phase ( F > fluid_isovalue)
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
DoubleArray CubeValues(2,2,2);
// int count_in=0,count_out=0;
// int nodx,nody,nodz;
// initialize lists for vertices for surfaces, common line
@ -1330,14 +1516,19 @@ int main(int argc, char **argv)
IntArray ws_tris(3,20);
// initialize list for line segments
IntArray nws_seg(2,20);
DTMutableList<Point> tmp(20);
DoubleArray Values(20);
DoubleArray ContactAngle(20);
DoubleArray Curvature(20);
DoubleArray InterfaceSpeed(20);
DoubleArray NormalVector(60);
// IntArray store;
int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=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)
// double s,s1,s2,s3; // Triangle sides (lengths)
Point A,B,C,P;
// double area;
@ -1507,6 +1698,24 @@ int main(int argc, char **argv)
dvc_UnpackValues(faceGrid, packThreads,dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, Phi, N);
//...................................................................................
//...........................................................................
// Copy the phase indicator field for the earlier timestep
dvc_Barrier();
dvc_CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
//...........................................................................
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
dvc_Barrier();
dvc_ComputePressure(nBlocks,nthreads,S,ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
//...........................................................................
int timestep = 0;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("No. of timesteps: %i \n", timestepMax);
@ -1515,7 +1724,6 @@ int main(int argc, char **argv)
// cudaStream_t stream;
// cudaStreamCreate(&stream);
//.......create and start timer............
double starttime,stoptime,cputime;
MPI_Barrier(MPI_COMM_WORLD);
@ -1523,9 +1731,16 @@ int main(int argc, char **argv)
//.........................................
sendtag = recvtag = 5;
if (rank==0) printf("-------------------------------------------------------------------\n");
if (rank==0) printf("timestep dEs sw Awn Ans Aws Lwns pw pn vwx vwy vwz vnx vny vnz \n");
if (rank==0) printf("-------------------------------------------------------------------\n");
if (rank==0){
printf("--------------------------------------------------------------------------------------\n");
printf("timestep dEs "); // Timestep, Change in Surface Energy
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");
}
//************ MAIN ITERATION LOOP ***************************************/
while (timestep < timestepMax){
@ -1770,7 +1985,6 @@ int main(int argc, char **argv)
// Wait for completion of D3Q7 communication
MPI_Waitall(6,req1,stat1);
MPI_Waitall(6,req2,stat2);
//...................................................................................
//...................................................................................
dvc_UnpackDenD3Q7(faceGrid,packThreads,dvcSendList_x,sendCount_x,sendbuf_x,2,Den,N);
@ -1780,7 +1994,6 @@ int main(int argc, char **argv)
dvc_UnpackDenD3Q7(faceGrid,packThreads,dvcSendList_Y,sendCount_Y,sendbuf_Y,2,Den,N);
dvc_UnpackDenD3Q7(faceGrid,packThreads,dvcSendList_Z,sendCount_Z,sendbuf_Z,2,Den,N);
//...................................................................................
//*************************************************************************
// Compute the phase indicator field and reset Copy, Den
@ -1886,8 +2099,16 @@ int main(int argc, char **argv)
// Iteration completed!
timestep++;
//...................................................................
if (timestep%1000 == 995){
//...........................................................................
// Copy the phase indicator field for the earlier timestep
dvc_Barrier();
dvc_CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
//...........................................................................
}
if (timestep%1000 == 0){
//...........................................................................
// Copy the data for for the analysis timestep
//...........................................................................
// Copy the phase from the GPU -> CPU
//...........................................................................
@ -1896,8 +2117,145 @@ int main(int argc, char **argv)
dvc_CopyToHost(Phase.data,Phi,N*sizeof(double));
dvc_CopyToHost(Press,Pressure,N*sizeof(double));
dvc_CopyToHost(Vel,Velocity,3*N*sizeof(double));
MPI_Barrier(MPI_COMM_WORLD);
}
if (timestep%1000 == 5){
//...........................................................................
// Copy the phase indicator field for the later timestep
dvc_Barrier();
dvc_CopyToHost(Phase_tminus.data,Phi,N*sizeof(double));
//...........................................................................
// Calculate the time derivative of the phase indicator field
for (n=0; n<N; n++) dPdt(n) = 0.1*(Phase_tplus(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);
//...........................................................................
// 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
//...........................................................................
// Mean Curvature
//...........................................................................
CommunicateMeshHalo(MeanCurvature, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
//...........................................................................
// Gaussian Curvature
//...........................................................................
CommunicateMeshHalo(GaussCurvature, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
// Gradient of the phase indicator field
//...........................................................................
CommunicateMeshHalo(Phase_x, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
CommunicateMeshHalo(Phase_y, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
CommunicateMeshHalo(Phase_z, MPI_COMM_WORLD,
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
sendMeshData_xZ,sendMeshData_Xz,sendMeshData_yz,sendMeshData_YZ,sendMeshData_yZ,sendMeshData_Yz,
recvMeshData_x,recvMeshData_y,recvMeshData_z,recvMeshData_X,recvMeshData_Y,recvMeshData_Z,
recvMeshData_xy,recvMeshData_XY,recvMeshData_xY,recvMeshData_Xy,recvMeshData_xz,recvMeshData_XZ,
recvMeshData_xZ,recvMeshData_Xz,recvMeshData_yz,recvMeshData_YZ,recvMeshData_yZ,recvMeshData_Yz,
sendList_x,sendList_y,sendList_z,sendList_X,sendList_Y,sendList_Z,
sendList_xy,sendList_XY,sendList_xY,sendList_Xy,sendList_xz,sendList_XZ,
sendList_xZ,sendList_Xz,sendList_yz,sendList_YZ,sendList_yZ,sendList_Yz,
sendCount_x,sendCount_y,sendCount_z,sendCount_X,sendCount_Y,sendCount_Z,
sendCount_xy,sendCount_XY,sendCount_xY,sendCount_Xy,sendCount_xz,sendCount_XZ,
sendCount_xZ,sendCount_Xz,sendCount_yz,sendCount_YZ,sendCount_yZ,sendCount_Yz,
recvList_x,recvList_y,recvList_z,recvList_X,recvList_Y,recvList_Z,
recvList_xy,recvList_XY,recvList_xY,recvList_Xy,recvList_xz,recvList_XZ,
recvList_xZ,recvList_Xz,recvList_yz,recvList_YZ,recvList_yZ,recvList_Yz,
recvCount_x,recvCount_y,recvCount_z,recvCount_X,recvCount_Y,recvCount_Z,
recvCount_xy,recvCount_XY,recvCount_xY,recvCount_Xy,recvCount_xz,recvCount_XZ,
recvCount_xZ,recvCount_Xz,recvCount_yz,recvCount_YZ,recvCount_yZ,recvCount_Yz,
rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z,rank_xy,rank_XY,rank_xY,
rank_Xy,rank_xz,rank_XZ,rank_xZ,rank_Xz,rank_yz,rank_YZ,rank_yZ,rank_Yz);
//...........................................................................
//...........................................................................
// Compute areas using porous medium marching cubes algorithm
// McClure, Adalsteinsson, et al. (2007)
@ -1906,17 +2264,17 @@ int main(int argc, char **argv)
nwp_volume = 0.0;
As = 0.0;
// Compute the normal vector
for (n=0; n<N; n++){
Norm_x[n] = 0.0;
Norm_y[n] = 0.0;
Norm_z[n] = 0.0;
}
// Compute phase averages
p_n = p_w = 0.0;
vx_w = vy_w = vz_w = 0.0;
vx_n = vy_n = vz_n = 0.0;
pan = paw = 0.0;
vaw(0) = vaw(1) = vaw(2) = 0.0;
van(0) = van(1) = van(2) = 0.0;
vawn(0) = vawn(1) = vawn(2) = 0.0;
Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
Gws(0) = Gws(1) = Gws(2) = 0.0;
Gws(3) = Gws(4) = Gws(5) = 0.0;
Gns(0) = Gns(1) = Gns(2) = 0.0;
Gns(3) = Gns(4) = Gns(5) = 0.0;
vol_w = vol_n =0.0;
for (c=0;c<ncubes;c++){
@ -1942,11 +2300,11 @@ int main(int argc, char **argv)
// volume the excludes the interfacial region
vol_n += 0.125;
// pressure
p_n += 0.125*Press[n];
pan += 0.125*Press[n];
// velocity
vx_n += 0.125*Vel[3*n];
vy_n += 0.125*Vel[3*n+1];
vz_n += 0.125*Vel[3*n+2];
van(0) += 0.125*Vel[3*n];
van(1) += 0.125*Vel[3*n+1];
van(2) += 0.125*Vel[3*n+2];
}
// volume averages over the wetting phase
@ -1954,230 +2312,48 @@ int main(int argc, char **argv)
// volume the excludes the interfacial region
vol_w += 0.125;
// pressure
p_w += 0.125*Press[n];
paw += 0.125*Press[n];
// velocity
vx_w += 0.125*Vel[3*n];
vy_w += 0.125*Vel[3*n+1];
vz_w += 0.125*Vel[3*n+2];
vaw(0) += 0.125*Vel[3*n];
vaw(1) += 0.125*Vel[3*n+1];
vaw(2) += 0.125*Vel[3*n+2];
}
}
}
//...........................................................................
// 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);
// 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;
//n_nw_tris_beg = 0;// n_nw_tris;
//n_ns_tris_beg = 0;//n_ns_tris;
//n_ws_tris_beg = 0;//n_ws_tris;
// Integrate the contact angle
efawns += pmmc_CubeContactAngle(CubeValues,Values,Phase_x,Phase_y,Phase_z,SignDist_x,SignDist_y,SignDist_z,
local_nws_pts,i,j,k,n_local_nws_pts);
// if there is a solid phase interface in the grid cell
if (Interface(SignDist,solid_isovalue,i,j,k) == 1){
/////////////////////////////////////////
/// CONSTRUCT THE LOCAL SOLID SURFACE ///
/////////////////////////////////////////
// find the local solid surface
// SOL_SURF(SignDist,0.0,Phase,fluid_isovalue,i,j,k, Nx,Ny,Nz,local_sol_pts,n_local_sol_pts,
// local_sol_tris,n_local_sol_tris,values);
// find the local solid surface using the regular Marching Cubes algorithm
SolidMarchingCubes(SignDist,0.0,Phase,fluid_isovalue,i,j,k,Nx,Ny,Nz,local_sol_pts,n_local_sol_pts,
local_sol_tris,n_local_sol_tris,values);
/////////////////////////////////////////
//////// TRIM THE SOLID SURFACE /////////
/////////////////////////////////////////
/* TRIM(local_sol_pts, n_local_sol_pts, fluid_isovalue,local_sol_tris, n_local_sol_tris,
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts,
Phase, SignDist, i, j, k, newton_steps);
*/
TRIM(local_sol_pts, n_local_sol_pts, fluid_isovalue,local_sol_tris, n_local_sol_tris,
ns_pts, n_ns_pts, ns_tris, n_ns_tris, ws_pts, n_ws_pts,
ws_tris, n_ws_tris, values, local_nws_pts, n_local_nws_pts);
/////////////////////////////////////////
//////// WRITE COMMON LINE POINTS ///////
//////// TO MAIN ARRAYS ///////
/////////////////////////////////////////
// SORT THE LOCAL COMMON LINE POINTS
/////////////////////////////////////////
// Make sure the first common line point is on a face
// Common curve points are located pairwise and must
// be searched and rearranged accordingly
for (p=0; p<n_local_nws_pts-1; p++){
P = local_nws_pts(p);
if ( P.x == 1.0*i || P.x ==1.0*(i+1)||
P.y == 1.0*j || P.y == 1.0*(j+1) ||
P.z == 1.0*k || P.z == 1.0*(k+1) ){
if (p%2 == 0){
// even points
// Swap the pair of points
local_nws_pts(p) = local_nws_pts(0);
local_nws_pts(0) = P;
P = local_nws_pts(p+1);
local_nws_pts(p+1) = local_nws_pts(1);
local_nws_pts(1) = P;
p = n_local_nws_pts;
}
else{
// odd points - flip the order
local_nws_pts(p) = local_nws_pts(p-1);
local_nws_pts(p-1) = P;
p-=2;
}
// guarantee exit from the loop
}
}
// Two common curve points per triangle
// 0-(1=2)-(3=4)-...
for (p=1; p<n_local_nws_pts-1; p+=2){
A = local_nws_pts(p);
for (q=p+1; q<n_local_nws_pts; q++){
B = local_nws_pts(q);
if ( A.x == B.x && A.y == B.y && A.z == B.z){
if (q%2 == 0){
// even points
// Swap the pair of points
local_nws_pts(q) = local_nws_pts(p+1);
local_nws_pts(p+1) = B;
B = local_nws_pts(q+1);
local_nws_pts(q+1) = local_nws_pts(p+2);
local_nws_pts(p+2) = B;
q = n_local_nws_pts;
}
else{
// odd points - flip the order
local_nws_pts(q) = local_nws_pts(q-1);
local_nws_pts(q-1) = B;
q-=2;
}
}
}
}
map = n_nws_pts = 0;
nws_pts(n_nws_pts++) = local_nws_pts(0);
for (p=2; p < n_local_nws_pts; p+=2){
nws_pts(n_nws_pts++) = local_nws_pts(p);
}
nws_pts(n_nws_pts++) = local_nws_pts(n_local_nws_pts-1);
for (q=0; q < n_nws_pts-1; q++){
nws_seg(0,n_nws_seg) = map+q;
nws_seg(1,n_nws_seg) = map+q+1;
n_nws_seg++;
}
// End of the common line sorting algorithm
/////////////////////////////////////////
/////////////////////////////////////////
////// CONSTRUCT THE nw SURFACE /////////
/////////////////////////////////////////
if ( n_local_nws_pts > 0){
n_nw_tris =0;
EDGE(Phase, fluid_isovalue, SignDist, i,j,k, Nx, Ny, Nz, nw_pts, n_nw_pts, nw_tris, n_nw_tris,
nws_pts, n_nws_pts);
}
else {
MC(Phase, fluid_isovalue, SignDist, i,j,k, nw_pts, n_nw_pts, nw_tris, n_nw_tris);
}
}
// Integrate the mean curvature
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
/////////////////////////////////////////
////// CONSTRUCT THE nw SURFACE /////////
/////////////////////////////////////////
else if (Fluid_Interface(Phase,SignDist,fluid_isovalue,i,j,k) == 1){
MC(Phase, fluid_isovalue, SignDist, i,j,k, nw_pts, n_nw_pts, nw_tris, n_nw_tris);
}
//******END OF BLOB PMMC*********************************************
//*******************************************************************
// Compute the Interfacial Areas, Common Line length for blob p
// nw surface
double temp;
for (r=0;r<n_nw_tris;r++){
A = nw_pts(nw_tris(0,r));
B = nw_pts(nw_tris(1,r));
C = nw_pts(nw_tris(2,r));
// 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));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) awn += sqrt(temp);
}
for (r=0;r<n_ns_tris;r++){
A = ns_pts(ns_tris(0,r));
B = ns_pts(ns_tris(1,r));
C = ns_pts(ns_tris(2,r));
// 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));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
//ans=ans+sqrt(s*(s-s1)*(s-s2)*(s-s3));
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) ans += sqrt(temp);
}
for (r=0;r<n_ws_tris;r++){
A = ws_pts(ws_tris(0,r));
B = ws_pts(ws_tris(1,r));
C = ws_pts(ws_tris(2,r));
// 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));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) aws += sqrt(temp);
}
for (r=0;r<n_local_sol_tris;r++){
A = local_sol_pts(local_sol_tris(0,r));
B = local_sol_pts(local_sol_tris(1,r));
C = local_sol_pts(local_sol_tris(2,r));
// 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));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0) As += sqrt(temp);
}
for (p=0; p < n_local_nws_pts-1; p++){
// Extract the line segment
A = local_nws_pts(p);
B = local_nws_pts(p+1);
// Compute the length of the segment
s = 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));
// Add the length to the common line
lwns += s;
}
//*******************************************************************
// Reset the triangle counts to zero
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;
//n_nw_tris_beg = 0;// n_nw_tris;
// n_ns_tris_beg = 0;//n_ns_tris;
// n_ws_tris_beg = 0;//n_ws_tris;
// n_nws_seg_beg = n_nws_seg;
//*******************************************************************
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);
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);
// 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);
lwns += pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
//...........................................................................
}
//...........................................................................
MPI_Barrier(MPI_COMM_WORLD);
@ -2187,17 +2363,19 @@ int main(int argc, char **argv)
MPI_Allreduce(&aws,&aws_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&lwns,&lwns_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&As,&As_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&Jwn,&Jwn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&efawns,&efawns_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
// Phase averages
MPI_Allreduce(&vol_w,&vol_w_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vol_n,&vol_n_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&p_w,&p_w_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&p_n,&p_n_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vx_w,&vx_w_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vy_w,&vy_w_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vz_w,&vz_w_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vx_n,&vx_n_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vy_n,&vy_n_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vz_n,&vz_n_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&paw,&paw_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&pan,&pan_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vaw(0),&vaw_global(0),3,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&van(0),&van_global(0),3,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&vawn(0),&vawn_global(0),3,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&Gwn(0),&Gwn_global(0),6,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&Gns(0),&Gns_global(0),6,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&Gws(0),&Gws_global(0),6,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
//.........................................................................
// Compute the change in the total surface energy based on the defined interval
@ -2208,18 +2386,28 @@ int main(int argc, char **argv)
dEs = 6.01603*alpha*(dAwn + 1.05332*Ps*dAns);
dAwn = -awn_global; // Get ready for the next analysis interval
dAns = -ans_global;
// Finalize the phase averages
// Normalize the phase averages
// (density of both components = 1.0)
p_w_global = p_w_global / vol_w_global;
vx_w_global = vx_w_global / vol_w_global;
vy_w_global = vy_w_global / vol_w_global;
vz_w_global = vz_w_global / vol_w_global;
p_n_global = p_n_global / vol_n_global;
vx_n_global = vx_n_global / vol_n_global;
vy_n_global = vy_n_global / vol_n_global;
vz_n_global = vz_n_global / vol_n_global;
paw_global = paw_global / vol_w_global;
vaw_global(0) = vaw_global(0) / vol_w_global;
vaw_global(1) = vaw_global(1) / vol_w_global;
vaw_global(2) = vaw_global(2) / vol_w_global;
pan_global = pan_global / vol_n_global;
van_global(0) = van_global(0) / vol_n_global;
van_global(1) = van_global(1) / vol_n_global;
van_global(2) = van_global(2) / vol_n_global;
// Normalize surface averages by the interfacial area
Jwn_global /= awn_global;
efawns_global /= lwns_global;
for (i=0; i<3; i++) vawn_global(i) /= awn_global;
for (i=0; i<6; i++) Gwn_global(i) /= awn_global;
for (i=0; i<6; i++) Gns_global(i) /= ans_global;
for (i=0; i<6; i++) Gws_global(i) /= aws_global;
sat_w = 1.0 - nwp_volume_global*iVol_global/porosity;
// Area and common curve terms
// Compute the specific interfacial areas and common line length (per unit volume)
awn_global = awn_global*iVol_global;
ans_global = ans_global*iVol_global;
aws_global = aws_global*iVol_global;
@ -2236,12 +2424,26 @@ int main(int argc, char **argv)
printf("Area ws = %f \n", aws_global);
printf("Change in surface energy = %f \n", dEs);
printf("-------------------------------- \n");
*/
printf("%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g \n",
timestep,dEs,sat_w,
awn_global,ans_global,aws_global, lwns_global, p_w_global, p_n_global,
vx_w_global, vy_w_global, vz_w_global,
vx_n_global, vy_n_global, vz_n_global);
*/
printf("%i %.5g ",timestep-5,dEs); // change in surface energy
printf("%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure
printf("%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase
printf("%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase
printf("%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas
printf("%.5g ",Jwn_global); // curvature of wn interface
printf("%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface
printf("%.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
printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
Gns_global(0),Gns_global(1),Gns_global(2),Gns_global(3),Gns_global(4),Gns_global(5)); // orientation of ns interface
printf("%.5g %.5g %.5g %.5g %.5g %.5g \n",
Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface
}
}
@ -2281,7 +2483,8 @@ int main(int argc, char **argv)
//#ifdef WriteOutput
FILE *PHASE;
PHASE = fopen(LocalRankFilename,"wb");
fwrite(Phase.data,8,N,PHASE);
// fwrite(Phase.data,8,N,PHASE);
fwrite(MeanCurvature.data,8,N,PHASE);
fclose(PHASE);
//#endif

12
include/Array.h
View File

@ -312,14 +312,14 @@ double DoubleArray::e(int i, int j, int k)
extern DoubleArray IncreaseSize(DoubleArray &A, int addLength)
{
if (addLength<0) {
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
printf("IncreaseSize(Array,Length) Length needs to be >0.");
return DoubleArray();
}
int newM,newN,newO;
if (A.o>1) {
if (addLength%(A.m*A.n)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m*n");
return DoubleArray();
}
newM = A.m;
@ -328,7 +328,7 @@ extern DoubleArray IncreaseSize(DoubleArray &A, int addLength)
}
else if (A.n>1) {
if (addLength%(A.m)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m");
return DoubleArray();
}
newM = A.m;
@ -348,14 +348,14 @@ extern DoubleArray IncreaseSize(DoubleArray &A, int addLength)
extern IntArray IncreaseSize(IntArray &A, int addLength)
{
if (addLength<0) {
printf("IncreaseSize(Array,Length)","Length needs to be >0.");
printf("IncreaseSize(Array,Length) Length needs to be >0.");
return IntArray();
}
int newM,newN,newO;
if (A.o>1) {
if (addLength%(A.m*A.n)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m*n");
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m*n");
return IntArray();
}
newM = A.m;
@ -364,7 +364,7 @@ extern IntArray IncreaseSize(IntArray &A, int addLength)
}
else if (A.n>1) {
if (addLength%(A.m)!=0) {
printf("IncreaseSize(Array,Length)","Length needs to be a multiple of m");
printf("IncreaseSize(Array,Length) Length needs to be a multiple of m");
return IntArray();
}
newM = A.m;

15
include/Domain.h
View File

@ -399,3 +399,18 @@ inline void ReadCheckpoint(char *FILENAME, double *cDen, double *cDistEven, doub
}
File.close();
}
inline void ReadBinaryFile(char *FILENAME, double *Data, int N)
{
int n;
double value;
ifstream File(FILENAME,ios::binary);
for (n=0; n<N; n++){
// Write the two density values
File.read((char*) &value, sizeof(value));
Data[n] = value;
}
File.close();
}

454
include/pmmc.h
View File

@ -3402,15 +3402,24 @@ inline void pmmc_MeshCurvature(DoubleArray &f, DoubleArray &MeanCurvature, Doubl
fyz = 0.25*(f(i,j+1,k+1) - f(i,j+1,k-1) - f(i,j-1,k+1) + f(i,j-1,k-1));
// Evaluate the Mean Curvature
denominator = pow(sqrt(fx*fx + fy*fy + fz*fz),3);
if (denominator == 0.0) denominator = 1.0;
MeanCurvature(i,j,k)=(1.0/denominator)*((fyy+fzz)*fx*fx + (fxx+fzz)*fy*fy + (fxx+fyy)*fz*fz
-2.0*fx*fy*fxy - 2.0*fx*fz*fxz - 2.0*fy*fz*fyz);
if (denominator == 0.0){
MeanCurvature(i,j,k) = 0.0;
}
else{
MeanCurvature(i,j,k)=(1.0/denominator)*((fyy+fzz)*fx*fx + (fxx+fzz)*fy*fy + (fxx+fyy)*fz*fz
-2.0*fx*fy*fxy - 2.0*fx*fz*fxz - 2.0*fy*fz*fyz);
}
// Evaluate the Gaussian Curvature
denominator = pow(fx*fx + fy*fy + fz*fz,2);
if (denominator == 0.0) denominator = 1.0;
GaussCurvature(i,j,k) = (1.0/denominator)*(fx*fx*(fyy*fzz-fyz*fyz) + fy*fy*(fxx*fzz-fxz*fxz) + fz*fz*(fxx*fyy-fxy*fxy)
+2.0*(fx*fy*(fxz*fyz-fxy*fzz) + fy*fz*(fxy*fxz-fyz*fxx)
+ fx*fz*(fxy*fyz-fxz*fyy)));
if (denominator == 0.0){
GaussCurvature(i,j,k) = 0.0;
}
else{
GaussCurvature(i,j,k) = (1.0/denominator)*(fx*fx*(fyy*fzz-fyz*fyz) + fy*fy*(fxx*fzz-fxz*fxz) + fz*fz*(fxx*fyy-fxy*fxy)
+2.0*(fx*fy*(fxz*fyz-fxy*fzz) + fy*fz*(fxy*fxz-fyz*fxx)
+ fx*fz*(fxy*fyz-fxz*fyy)));
}
}
}
}
@ -3440,6 +3449,7 @@ inline void pmmc_CubeListFromBlobs()
}
//--------------------------------------------------------------------------------------------------------
inline double pmmc_CubeSurfaceArea(DTMutableList<Point> &Points, IntArray &Triangles, int ntris)
{
int r;
@ -3461,7 +3471,25 @@ inline double pmmc_CubeSurfaceArea(DTMutableList<Point> &Points, IntArray &Trian
return area;
}
//--------------------------------------------------------------------------------------------------------
inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
inline double pmmc_CubeCurveLength(DTMutableList<Point> &Points, int npts)
{
int p;
double s,lwns;
Point A,B;
lwns = 0.0;
for (p=0; p < npts-1; p++){
// Extract the line segment
A = Points(p);
B = Points(p+1);
// Compute the length of the segment
s = 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));
// Add the length to the common line
lwns += s;
}
return lwns;
}
inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DoubleArray &MeshValues, DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &SurfaceValues, int i, int j, int k, int npts, int ntris)
{
Point A,B,C;
@ -3472,6 +3500,16 @@ inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList
double a,b,c,d,e,f,g,h;
double integral;
// Copy the curvature values for the cube
CubeValues(0,0,0) = MeshValues(i,j,k);
CubeValues(1,0,0) = MeshValues(i+1,j,k);
CubeValues(0,1,0) = MeshValues(i,j+1,k);
CubeValues(1,1,0) = MeshValues(i+1,j+1,k);
CubeValues(0,0,1) = MeshValues(i,j,k+1);
CubeValues(1,0,1) = MeshValues(i+1,j,k+1);
CubeValues(0,1,1) = MeshValues(i,j+1,k+1);
CubeValues(1,1,1) = MeshValues(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
// Evaluate the coefficients
a = CubeValues(0,0,0);
@ -3511,13 +3549,13 @@ inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList
}
//--------------------------------------------------------------------------------------------------------
inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &CurveValues,
DTMutableList<Point> &Points, int npts)
DTMutableList<Point> &Points, int i, int j, int k, int npts)
{
int p;
Point A,B;
double vA,vB;
double x,y,z;
double s,s1,s2,s3,temp;
// double s,s1,s2,s3,temp;
double a,b,c,d,e,f,g,h;
double integral;
double length;
@ -3535,7 +3573,10 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
for (p=0; p<npts; p++){
A = Points(p);
CurveValues(p) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
integral = 0.0;
@ -3543,9 +3584,11 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
// Extract the line segment
A = Points(p);
B = Points(p+1);
vA = CurveValues(p);
vB = CurveValues(p+1);
// Compute the length of the segment
length = 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));
integral += 0.5*length*(CurveValues(p) + CurveValues(p+1));
integral += 0.5*length*(vA + vB);
}
return integral;
}
@ -3558,7 +3601,8 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
int p;
Point A,B;
double vA,vB;
double s,s1,s2,s3,temp;
double x,y,z;
// double s,s1,s2,s3,temp;
double a,b,c,d,e,f,g,h;
double integral;
double length;
@ -3600,9 +3644,12 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (int p=0; p<npts; p++){
for (p=0; p<npts; p++){
A = Points(p);
CurveValues(p) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
integral = 0.0;
@ -3610,10 +3657,381 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
// Extract the line segment
A = Points(p);
B = Points(p+1);
vA = CurveValues(p);
vB = CurveValues(p+1);
// Compute the length of the segment
length = 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));
integral += 0.5*length*(CurveValues(p) + CurveValues(p+1));
integral += 0.5*length*(vA + vB);
}
return integral;
}//--------------------------------------------------------------------------------------------------------
}
//--------------------------------------------------------------------------------------------------------
inline void pmmc_CubeSurfaceInterpVector(DoubleArray &Vec_x, DoubleArray &Vec_y, DoubleArray &Vec_z,
DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &SurfaceVector, DoubleArray &VecAvg, int i, int j, int k, int npts, int ntris)
{
Point A,B,C;
int p;
double vA,vB,vC;
double x,y,z;
double s,s1,s2,s3,temp;
double a,b,c,d,e,f,g,h;
// ................x component .............................
// Copy the curvature values for the cube
CubeValues(0,0,0) = Vec_x(i,j,k);
CubeValues(1,0,0) = Vec_x(i+1,j,k);
CubeValues(0,1,0) = Vec_x(i,j+1,k);
CubeValues(1,1,0) = Vec_x(i+1,j+1,k);
CubeValues(0,0,1) = Vec_x(i,j,k+1);
CubeValues(1,0,1) = Vec_x(i+1,j,k+1);
CubeValues(0,1,1) = Vec_x(i,j+1,k+1);
CubeValues(1,1,1) = Vec_x(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
// ................y component .............................
// Copy the curvature values for the cube
CubeValues(0,0,0) = Vec_y(i,j,k);
CubeValues(1,0,0) = Vec_y(i+1,j,k);
CubeValues(0,1,0) = Vec_y(i,j+1,k);
CubeValues(1,1,0) = Vec_y(i+1,j+1,k);
CubeValues(0,0,1) = Vec_y(i,j,k+1);
CubeValues(1,0,1) = Vec_y(i+1,j,k+1);
CubeValues(0,1,1) = Vec_y(i,j+1,k+1);
CubeValues(1,1,1) = Vec_y(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(npts+p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
// ................z component .............................
// Copy the curvature values for the cube
CubeValues(0,0,0) = Vec_z(i,j,k);
CubeValues(1,0,0) = Vec_z(i+1,j,k);
CubeValues(0,1,0) = Vec_z(i,j+1,k);
CubeValues(1,1,0) = Vec_z(i+1,j+1,k);
CubeValues(0,0,1) = Vec_z(i,j,k+1);
CubeValues(1,0,1) = Vec_z(i+1,j,k+1);
CubeValues(0,1,1) = Vec_z(i,j+1,k+1);
CubeValues(1,1,1) = Vec_z(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(2*npts+p) = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
//.............................................................................
for (int r=0; r<ntris; r++){
A = Points(Triangles(0,r));
B = Points(Triangles(1,r));
C = Points(Triangles(2,r));
// 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));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0){
// Increment the averaged values
// x component
vA = SurfaceVector(Triangles(0,r));
vB = SurfaceVector(Triangles(1,r));
vC = SurfaceVector(Triangles(2,r));
VecAvg(0) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
// y component
vA = SurfaceVector(npts+Triangles(0,r));
vB = SurfaceVector(npts+Triangles(1,r));
vC = SurfaceVector(npts+Triangles(2,r));
VecAvg(1) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
// z component
vA = SurfaceVector(2*npts+Triangles(0,r));
vB = SurfaceVector(2*npts+Triangles(1,r));
vC = SurfaceVector(2*npts+Triangles(2,r));
VecAvg(2) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
}
}
}
//--------------------------------------------------------------------------------------------------------
inline double pmmc_CubeSurfaceOrientation(DoubleArray &Orientation, DTMutableList<Point> &Points, IntArray &Triangles, int ntris)
{
int r;
double temp,area,s,s1,s2,s3;
double nx,ny,nz,normsq;
Point A,B,C;
area = 0.0;
for (r=0;r<ntris;r++){
A = Points(Triangles(0,r));
B = Points(Triangles(1,r));
C = Points(Triangles(2,r));
// Compute the triangle normal vector
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);
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));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0){
temp = sqrt(temp);
area += temp;
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;
}
//--------------------------------------------------------------------------------------------------------
inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray &P_y, DoubleArray &P_z,
DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
DoubleArray &SurfaceVector, DoubleArray &SurfaceValues, DoubleArray &AvgVel,
int i, int j, int k, int npts, int ntris)
{
Point A,B,C;
int p;
double vA,vB,vC;
double vAx,vBx,vCx,vAy,vBy,vCy,vAz,vBz,vCz;
double x,y,z;
double s,s1,s2,s3,temp;
double a,b,c,d,e,f,g,h;
double norm, zeta;
// ................x component .............................
// Copy the x derivative values for the cube
CubeValues(0,0,0) = P_x(i,j,k);
CubeValues(1,0,0) = P_x(i+1,j,k);
CubeValues(0,1,0) = P_x(i,j+1,k);
CubeValues(1,1,0) = P_x(i+1,j+1,k);
CubeValues(0,0,1) = P_x(i,j,k+1);
CubeValues(1,0,1) = P_x(i+1,j,k+1);
CubeValues(0,1,1) = P_x(i,j+1,k+1);
CubeValues(1,1,1) = P_x(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
// ................y component .............................
// Copy the y derivative values for the cube
CubeValues(0,0,0) = P_y(i,j,k);
CubeValues(1,0,0) = P_y(i+1,j,k);
CubeValues(0,1,0) = P_y(i,j+1,k);
CubeValues(1,1,0) = P_y(i+1,j+1,k);
CubeValues(0,0,1) = P_y(i,j,k+1);
CubeValues(1,0,1) = P_y(i+1,j,k+1);
CubeValues(0,1,1) = P_y(i,j+1,k+1);
CubeValues(1,1,1) = P_y(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(npts+p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
// ................z component .............................
// Copy the z derivative values for the cube
CubeValues(0,0,0) = P_z(i,j,k);
CubeValues(1,0,0) = P_z(i+1,j,k);
CubeValues(0,1,0) = P_z(i,j+1,k);
CubeValues(1,1,0) = P_z(i+1,j+1,k);
CubeValues(0,0,1) = P_z(i,j,k+1);
CubeValues(1,0,1) = P_z(i+1,j,k+1);
CubeValues(0,1,1) = P_z(i,j+1,k+1);
CubeValues(1,1,1) = P_z(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
SurfaceVector(2*npts+p) = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
}
//.............................................................................
// Compute the normal and the speed at points on the interface
// Copy the curvature values for the cube
CubeValues(0,0,0) = dPdt(i,j,k);
CubeValues(1,0,0) = dPdt(i+1,j,k);
CubeValues(0,1,0) = dPdt(i,j+1,k);
CubeValues(1,1,0) = dPdt(i+1,j+1,k);
CubeValues(0,0,1) = dPdt(i,j,k+1);
CubeValues(1,0,1) = dPdt(i+1,j,k+1);
CubeValues(0,1,1) = dPdt(i,j+1,k+1);
CubeValues(1,1,1) = dPdt(i+1,j+1,k+1);
// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
a = CubeValues(0,0,0);
b = CubeValues(1,0,0)-a;
c = CubeValues(0,1,0)-a;
d = CubeValues(0,0,1)-a;
e = CubeValues(1,1,0)-a-b-c;
f = CubeValues(1,0,1)-a-b-d;
g = CubeValues(0,1,1)-a-c-d;
h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
for (p=0; p<npts; p++){
A = Points(p);
// evaluate time derivative on the surface
x = A.x-1.0*i;
y = A.y-1.0*j;
z = A.z-1.0*k;
zeta = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
// compute the normal
x = SurfaceVector(p);
y = SurfaceVector(npts+p);
z = SurfaceVector(2*npts+p);
norm = sqrt(x*x+y*y+z*z);
// Save the surface values and normal vector
if (norm > 0.0){
SurfaceValues(p) = -zeta/norm;
SurfaceVector(p) = x/norm;
SurfaceVector(npts+p) = y/norm;
SurfaceVector(2*npts+p) = z/norm;
}
else{
SurfaceValues(p) = 0.0;
SurfaceVector(p) = 0.0;
SurfaceVector(npts+p) = 0.0;
SurfaceVector(2*npts+p) = 0.0;
}
}
//.............................................................................
// Compute the average speed of the interface
for (int r=0; r<ntris; r++){
A = Points(Triangles(0,r));
B = Points(Triangles(1,r));
C = Points(Triangles(2,r));
// 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));
s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
s = 0.5*(s1+s2+s3);
temp = s*(s-s1)*(s-s2)*(s-s3);
if (temp > 0.0){
// Surface value (speed)
vA = SurfaceValues(Triangles(0,r));
vB = SurfaceValues(Triangles(1,r));
vC = SurfaceValues(Triangles(2,r));
// Increment the averaged values
// x component
vAx = SurfaceVector(Triangles(0,r))*vA;
vBx = SurfaceVector(Triangles(1,r))*vB;
vCx = SurfaceVector(Triangles(2,r))*vC;
// y component
vAy = SurfaceVector(npts+Triangles(0,r))*vA;
vBy = SurfaceVector(npts+Triangles(1,r))*vB;
vCy = SurfaceVector(npts+Triangles(2,r))*vC;
// z component
vAz = SurfaceVector(2*npts+Triangles(0,r))*vA;
vBz = SurfaceVector(2*npts+Triangles(1,r))*vB;
vCz = SurfaceVector(2*npts+Triangles(2,r))*vC;
AvgVel(0) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx);
AvgVel(1) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy);
AvgVel(2) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz);
// Update the Averages. Differentiate between advancing (0,1,2) and receding (3,4,5) interfaces
// All points on a triangle have the same orientation in the color gradient
/* if (vA > 0.0){
// Advancing interface
AvgVel(0) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx);
AvgVel(1) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy);
AvgVel(2) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz);
}
else{
// Receding interface
AvgVel(3) += sqrt(temp)*0.33333333333333333*(vAx+vBx+vCx);
AvgVel(4) += sqrt(temp)*0.33333333333333333*(vAy+vBy+vCy);
AvgVel(5) += sqrt(temp)*0.33333333333333333*(vAz+vBz+vCz);
}
*/
}
}
//.............................................................................
}
//--------------------------------------------------------------------------------------------------------

436
pmmc/Analysis.cpp Normal file
View File

@ -0,0 +1,436 @@
#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,N;
int i,j,k,p,q,r,n;
int nspheres;
double Lx,Ly,Lz;
//.......................................................................
Nx = Ny = Nz = 60;
cout << "Enter Domain size " << endl;
cout << "Nx = " << endl;
cin >> Nx;
Ny = Nz = Nx;
N = Nx*Ny*Nz;
//.......................................................................
// 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);
double iVol = 1.0/Nx/Ny/Nz;
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;
}
}
}
awn = aws = ans = lwns = 0.0;
nwp_volume = 0.0;
As = 0.0;
Jwn = 0.0;
efawns = 0.0;
// Compute phase averages
pan = paw = 0.0;
vaw(0) = vaw(1) = vaw(2) = 0.0;
van(0) = van(1) = van(2) = 0.0;
vawn(0) = vawn(1) = vawn(2) = 0.0;
Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
Gws(0) = Gws(1) = Gws(2) = 0.0;
Gws(3) = Gws(4) = Gws(5) = 0.0;
Gns(0) = Gns(1) = Gns(2) = 0.0;
Gns(3) = Gns(4) = Gns(5) = 0.0;
vol_w = vol_n =0.0;
// Read the input files for the phase, distance and pressure field
char PHASEFILE[16];
sprintf(PHASEFILE,"Phase.in");
ReadBinaryFile(PHASEFILE,Phase.data,Nx*Ny*Nz);
char DISTFILE[16];
sprintf(DISTFILE,"SignDist.in");
ReadBinaryFile(DISTFILE,SignDist.data,Nx*Ny*Nz);
/* FILE *PRESS
PRESS = fopen("Pressure.in","wb");
fread(Phase.data,8,N,PRESS);
fclose(PRESS);
FILE *VEL;
VEL = fopen("Pressure.in","wb");
fread(Phase.data,8,3*N,VEL);
fclose(VEL);
*/ //...........................................................................
// 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;
FILE *WN_TRIS;
WN_TRIS = fopen("wn-tris.out","w");
FILE *NS_TRIS;
NS_TRIS = fopen("ns-tris.out","w");
FILE *WS_TRIS;
WS_TRIS = fopen("ws-tris.out","w");
FILE *WNS_PTS;
WNS_PTS = fopen("wns-pts.out","w");
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);
// 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);
//.......................................................................................
// Write the triangle lists to text file
for (r=0;r<n_nw_tris;r++){
A = nw_pts(nw_tris(0,r));
B = nw_pts(nw_tris(1,r));
C = nw_pts(nw_tris(2,r));
fprintf(WN_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
}
for (r=0;r<n_ws_tris;r++){
A = ws_pts(ws_tris(0,r));
B = ws_pts(ws_tris(1,r));
C = ws_pts(ws_tris(2,r));
fprintf(WS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
}
for (r=0;r<n_ns_tris;r++){
A = ns_pts(ns_tris(0,r));
B = ns_pts(ns_tris(1,r));
C = ns_pts(ns_tris(2,r));
fprintf(NS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
}
for (p=0; p < n_nws_pts; p++){
P = nws_pts(p);
fprintf(WNS_PTS,"%f %f %f \n",P.x, P.y, P.z);
}
}
fclose(WN_TRIS);
fclose(NS_TRIS);
fclose(WS_TRIS);
fclose(WNS_PTS);
printf("Jwn = %f \n",Jwn);
printf("awn = %f \n",awn);
printf("efawns = %f \n",efawns);
printf("lwns = %f \n",lwns);
printf("efawns = %f \n",efawns/lwns);
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;
awn = awn*iVol;
aws = aws*iVol;
ans = ans*iVol;
lwns = lwns*iVol;
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 ", lwns); // common curve length
printf("%.5g ",efawns); // average contact angle
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 \n",
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.out","wb");
fwrite(Phase.data,8,N,PHASE);
fclose(PHASE);
FILE *SOLID;
SOLID = fopen("Distance.out","wb");
fwrite(SignDist.data,8,N,SOLID);
fclose(SOLID);
}

179
tests/TestContactAngle.cpp Normal file
View File

@ -0,0 +1,179 @@
#include <iostream>
#include <math.h>
#include "pmmc.h"
//#include "PointList.h"
//#include "Array.h"
#define RADIUS 15
#define CAPRAD 20
#define HEIGHT 15.5
#define N 60
#define PI 3.14159
int main (int argc, char *argv[])
{
// printf("Radius = %s \n,"RADIUS);
int SIZE = N*N*N;
int Nx,Ny,Nz;
Nx = Ny = Nz = N;
int i,j,k,p,q,r;
// double *Solid; // cylinder
// double *Phase; // region of the cylinder
// Solid = new double [SIZE];
// Phase = new double [SIZE];
DoubleArray SignDist(Nx,Ny,Nz);
DoubleArray Phase(Nx,Ny,Nz);
DoubleArray Fx(Nx,Ny,Nz);
DoubleArray Fy(Nx,Ny,Nz);
DoubleArray Fz(Nx,Ny,Nz);
DoubleArray Sx(Nx,Ny,Nz);
DoubleArray Sy(Nx,Ny,Nz);
DoubleArray Sz(Nx,Ny,Nz);
double fluid_isovalue = 0.0;
double solid_isovalue = 0.0;
/* ****************************************************************
VARIABLES FOR THE PMMC ALGORITHM
****************************************************************** */
//...........................................................................
// Averaging variables
//...........................................................................
double awn,ans,aws,lwns,nwp_volume;
double efawns;
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 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);
// IntArray store;
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;
// 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);
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;
Cx = Cy = Cz = N*0.51;
for (k=0; k<N; k++){
for (j=0; j<N; j++){
for (i=0; i<N; 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;
}
}
}
pmmc_MeshGradient(Phase,Fx,Fy,Fz,Nx,Ny,Nz);
pmmc_MeshGradient(SignDist,Sx,Sy,Sz,Nx,Ny,Nz);
// End of the loop to set the values
awn = aws = ans = lwns = 0.0;
nwp_volume = 0.0;
As = 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);
for (p=0;p<8;p++){
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0
&& SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
nwp_volume += 0.125;
}
}
// 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);
efawns += pmmc_CubeContactAngle(CubeValues,ContactAngle,Fx,Fy,Fz,Sx,Sy,Sz,local_nws_pts,i,j,k,n_local_nws_pts);
//*******************************************************************
// 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);
}
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("-------------------------------- \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);
*/
}

1
tests/TestCylinderAreas.cpp
View File

@ -145,7 +145,6 @@ int main (int argc, char *argv[])
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 Interfacial Areas, Common Line length for blob p

241
tests/TestInterfaceSpeed.cpp Normal file
View File

@ -0,0 +1,241 @@
#include <iostream>
#include <math.h>
#include "pmmc.h"
//#include "PointList.h"
//#include "Array.h"
#define RADIUS 15
#define CAPRAD 20
#define HEIGHT 15.5
#define N 60
#define SPEED -1
#define PI 3.14159
int main (int argc, char *argv[])
{
// printf("Radius = %s \n,"RADIUS);
int SIZE = N*N*N;
int Nx,Ny,Nz;
Nx = Ny = Nz = N;
int i,j,k,p,q,r;
// double *Solid; // cylinder
// double *Phase; // region of the cylinder
// Solid = new double [SIZE];
// Phase = new double [SIZE];
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 GaussCurvature(Nx,Ny,Nz);
DoubleArray MeanCurvature(Nx,Ny,Nz);
double fluid_isovalue = 0.0;
double solid_isovalue = 0.0;
/* ****************************************************************
VARIABLES FOR THE PMMC ALGORITHM
****************************************************************** */
//...........................................................................
// Averaging variables
//...........................................................................
double awn,ans,aws,lwns,nwp_volume;
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 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);
// IntArray store;
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;
// 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 wn_curvature(20);
DoubleArray InterfaceSpeed(20);
DoubleArray NormalVector(60);
DoubleArray vawn(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 = N*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;
dist2 = fabs(Cz-k)-HEIGHT;
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;
dist2 = fabs(Cz-k)-HEIGHT;
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;
dist2 = fabs(Cz-k)-HEIGHT;
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);
double norm;
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
norm = Phase_x(i,j,k)*Phase_x(i,j,k)+Phase_y(i,j,k)*Phase_y(i,j,k)+Phase_z(i,j,k)*Phase_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;
for (c=0;c<ncubes;c++){
// Get cube from the list
i = cubeList(0,c);
j = cubeList(1,c);
k = cubeList(2,c);
for (p=0;p<8;p++){
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0
&& SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
nwp_volume += 0.125;
}
}
// 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);
efawns += pmmc_CubeContactAngle(CubeValues,ContactAngle,Phase_x,Phase_y,Phase_z,Sx,Sy,Sz,local_nws_pts,i,j,k,n_local_nws_pts);
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature, nw_pts, nw_tris,
wn_curvature, 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);
// if (n_nw_pts>0) printf("speed %f \n",InterfaceSpeed(0));
//*******************************************************************
// 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<6;i++) vawn(i) /= awn;
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("Advancing Interface Velocity = %f,%f,%f \n",vawn(0),vawn(1),vawn(2));
printf("Receding Interface Velocity = %f,%f,%f \n",vawn(3),vawn(4),vawn(5));
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);
*/
}

16
tests/TestSphereCurvature.cpp
View File

@ -121,18 +121,8 @@ int main(int argc, char **argv)
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);
// Copy the curvature values for the cube
CubeValues(0,0,0) = MeanCurvature(i,j,k);
CubeValues(1,0,0) = MeanCurvature(i+1,j,k);
CubeValues(0,1,0) = MeanCurvature(i,j+1,k);
CubeValues(1,1,0) = MeanCurvature(i+1,j+1,k);
CubeValues(0,0,1) = MeanCurvature(i,j,k+1);
CubeValues(1,0,1) = MeanCurvature(i+1,j,k+1);
CubeValues(0,1,1) = MeanCurvature(i,j+1,k+1);
CubeValues(1,1,1) = MeanCurvature(i+1,j+1,k+1);
// Interpolate the curvature onto the surface
wn_curvature_sum += pmmc_CubeSurfaceInterpValue(CubeValues, nw_pts, nw_tris,
wn_curvature_sum += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature, nw_pts, nw_tris,
wn_curvature, i, j, k, n_nw_pts, n_nw_tris);
wn_area_sum += pmmc_CubeSurfaceArea(nw_pts, nw_tris, n_nw_tris);
@ -141,7 +131,7 @@ int main(int argc, char **argv)
printf("Mean Curvature Average = %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
FILE *CURVATURE;
/* FILE *CURVATURE;
CURVATURE = fopen("Curvature.dat","wb");
fwrite(MeanCurvature.data,8,N,CURVATURE);
fclose(CURVATURE);
@ -150,5 +140,5 @@ int main(int argc, char **argv)
DISTANCE = fopen("SignDist.dat","wb");
fwrite(Phase.data,8,N,DISTANCE);
fclose(DISTANCE);
*/
}

140
tests/TestSphereOrientation.cpp Normal file
View File

@ -0,0 +1,140 @@
#include <iostream>
#include <math.h>
#include "pmmc.h"
#include "Domain.h"
using namespace std;
/*
* Compare the measured and analytical curvature for a sphere
*
*/
int main(int argc, char **argv)
{
int i,j,k;
int Nx,Ny,Nz,N;
double Lx,Ly,Lz;
double fluid_isovalue=0.0;
double solid_isovalue=0.0;
Lx = Ly = Lz = 1.0;
Nx = Ny = Nz = 102;
N = Nx*Ny*Nz;
//...........................................................................
// Set up the cube list
//...........................................................................
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
IntArray cubeList(3,ncubes);
pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
//...........................................................................
//****************************************************************************************
// Create the structures needed to carry out the PMMC algorithm
//****************************************************************************************
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);
DoubleArray wn_curvature(20);
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;
// 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;
//****************************************************************************************
int nspheres=1;
// Set up the signed distance function for a sphere
double *cx,*cy,*cz,*rad;
cx = new double[nspheres];
cy = new double[nspheres];
cz = new double[nspheres];
rad = new double[nspheres];
// Sphere in the middle of the domain
cx[0] = cy[0] = cz[0] = 0.5;
rad[0] = 0.3;
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 CubeValues(2,2,2);
DoubleArray Orientation(6);
// Compute the signed distance function
SignedDistance(Phase.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
for (k=0; k<Nz; k++){
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
SignDist(i,j,k) = 100.0;
Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*j-0.5*Ny)*(1.0*j-0.5*Ny)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.3*Nx;
}
}
}
SignedDistance(SignDist.data,0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
double Gxx_sum,Gyy_sum,Gzz_sum,Gxy_sum,Gxz_sum,Gyz_sum;
double wn_curvature_sum = 0.0;
double wn_area_sum = 0.0;
Orientation(0) = Orientation(1) = Orientation(2) = 0.0;
Orientation(3) = Orientation(4) = Orientation(5) = 0.0;
pmmc_MeshGradient(Phase, Phase_x, Phase_y, Phase_z, Nx, Ny, Nz);
for (int 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);
wn_area_sum += pmmc_CubeSurfaceOrientation(Orientation,nw_pts,nw_tris,n_nw_tris);
}
printf("Gxx = %f, Analytical = 1/3 \n", Orientation(0)/wn_area_sum);
printf("Gyy = %f, Analytical = 1/3 \n", Orientation(1)/wn_area_sum);
printf("Gzz = %f, Analytical = 1/3 \n", Orientation(2)/wn_area_sum);
printf("Gxy = %f, Analytical = 0 \n", Orientation(3)/wn_area_sum);
printf("Gxz = %f, Analytical = 0 \n", Orientation(4)/wn_area_sum);
printf("Gyz = %f, Analytical = 0 \n", Orientation(5)/wn_area_sum);
}