Updated tests/lbpm_color_simulator to use averaging object
This commit is contained in:
James E McClure
parent
434bdf374e
commit
755cdecf2e
@ -360,7 +360,6 @@ void TwoPhase::SetupCubes(Domain &Dm){
|
||||
kfinish = Nz-1;
|
||||
if (Dm.BoundaryCondition !=0 && Dm.kproc==0) kstart = 4;
|
||||
if (Dm.BoundaryCondition !=0 && Dm.kproc==Dm.nprocz-1) kfinish = Nz-4;
|
||||
|
||||
nc=0;
|
||||
for (k=kstart; k<kfinish; k++){
|
||||
for (j=1; j<Ny-1; j++){
|
||||
@ -450,7 +449,6 @@ void TwoPhase::UpdateSolid(){
|
||||
}
|
||||
|
||||
void TwoPhase::UpdateMeshValues(){
|
||||
|
||||
//...........................................................................
|
||||
// Compute the gradients of the phase indicator and signed distance fields
|
||||
pmmc_MeshGradient(SDs,SDs_x,SDs_y,SDs_z,Nx,Ny,Nz);
|
||||
@ -461,10 +459,8 @@ void TwoPhase::UpdateMeshValues(){
|
||||
//...........................................................................
|
||||
// Update the time derivative of non-dimensional density field
|
||||
// Map Phase_tplus and Phase_tminus
|
||||
|
||||
for (int n=0; n<Nx*Ny*Nz; n++) dPdt(n) = 0.1*(Phase_tplus(n) - Phase_tminus(n));
|
||||
//...........................................................................
|
||||
//...........................................................................
|
||||
CommunicateMeshHalo(Vel_x, Dm.Comm,
|
||||
sendMeshData_x,sendMeshData_y,sendMeshData_z,sendMeshData_X,sendMeshData_Y,sendMeshData_Z,
|
||||
sendMeshData_xy,sendMeshData_XY,sendMeshData_xY,sendMeshData_Xy,sendMeshData_xz,sendMeshData_XZ,
|
||||
|
||||
@ -1050,13 +1050,6 @@ int main(int argc, char **argv)
|
||||
// Copy signed distance for device initialization
|
||||
CopyToDevice(dvcSDs, Averages.SDs.data, dist_mem_size);
|
||||
//...........................................................................
|
||||
// Phase indicator (in array form as needed by PMMC algorithm)
|
||||
DoubleArray Phase(Nx,Ny,Nz);
|
||||
|
||||
// Extra copies of phi needed to compute time derivatives on CPU
|
||||
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;
|
||||
@ -1230,7 +1223,7 @@ int main(int argc, char **argv)
|
||||
//...........................................................................
|
||||
// Copy the phase indicator field for the earlier timestep
|
||||
DeviceBarrier();
|
||||
CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
|
||||
CopyToHost(Averages.Phase_tplus.data,Phi,N*sizeof(double));
|
||||
//...........................................................................
|
||||
//...........................................................................
|
||||
// Copy the data for for the analysis timestep
|
||||
|
||||
@ -6,16 +6,15 @@
|
||||
#include <stdexcept>
|
||||
#include <fstream>
|
||||
|
||||
#include "pmmc.h"
|
||||
#include "Domain.h"
|
||||
#include "Extras.h"
|
||||
#include "D3Q19.h"
|
||||
#include "D3Q7.h"
|
||||
#include "Color.h"
|
||||
#include "Communication.h"
|
||||
#include "TwoPhase.h"
|
||||
#include "common/MPI.h"
|
||||
|
||||
|
||||
//#define WRITE_SURFACES
|
||||
|
||||
/*
|
||||
@ -276,7 +275,13 @@ int main(int argc, char **argv)
|
||||
printf("********************************************************\n");
|
||||
}
|
||||
|
||||
InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc,
|
||||
// Initialized domain and averaging framework for Two-Phase Flow
|
||||
int BC=pBC;
|
||||
Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
|
||||
Dm.CommInit(MPI_COMM_WORLD);
|
||||
TwoPhase Averages(Dm);
|
||||
|
||||
InitializeRanks( rank, nprocx, nprocy, nprocz, iproc, jproc, kproc,
|
||||
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 );
|
||||
@ -314,7 +319,6 @@ int main(int argc, char **argv)
|
||||
if (pBC) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
|
||||
double porosity, pore_vol;
|
||||
//...........................................................................
|
||||
DoubleArray SignDist(Nx,Ny,Nz);
|
||||
if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
|
||||
|
||||
//.......................................................................
|
||||
@ -322,7 +326,7 @@ int main(int argc, char **argv)
|
||||
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
|
||||
// WriteLocalSolidID(LocalRankFilename, id, N);
|
||||
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
|
||||
ReadBinaryFile(LocalRankFilename, SignDist.data, N);
|
||||
ReadBinaryFile(LocalRankFilename, Averages.SDs.data, N);
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
if (rank == 0) cout << "Domain set." << endl;
|
||||
|
||||
@ -343,11 +347,11 @@ int main(int argc, char **argv)
|
||||
for ( j=1;j<Ny-1;j++){
|
||||
for ( i=1;i<Nx-1;i++){
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
if (SignDist.data[n] > 0.0){
|
||||
if (Averages.SDs.data[n] > 0.0){
|
||||
id[n] = 2;
|
||||
}
|
||||
// compute the porosity (actual interface location used)
|
||||
if (SignDist.data[n] > 0.0){
|
||||
if (Averages.SDs.data[n] > 0.0){
|
||||
sum++;
|
||||
}
|
||||
}
|
||||
@ -389,7 +393,7 @@ int main(int argc, char **argv)
|
||||
for (i=0;i<Nx;i++){
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
id[n] = 1;
|
||||
SignDist.data[n] = max(SignDist.data[n],1.0*(2.5-k));
|
||||
Averages.SDs.data[n] = max(Averages.SDs.data[n],1.0*(2.5-k));
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -400,7 +404,7 @@ int main(int argc, char **argv)
|
||||
for (i=0;i<Nx;i++){
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
id[n] = 2;
|
||||
SignDist.data[n] = max(SignDist.data[n],1.0*(k-Nz+2.5));
|
||||
Averages.SDs.data[n] = max(Averages.SDs.data[n],1.0*(k-Nz+2.5));
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -958,146 +962,17 @@ int main(int argc, char **argv)
|
||||
AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
|
||||
AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
|
||||
AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
|
||||
// Copy signed distance for device initialization
|
||||
CopyToDevice(dvcSignDist, SignDist.data, dist_mem_size);
|
||||
//...........................................................................
|
||||
// Phase indicator (in array form as needed by PMMC algorithm)
|
||||
DoubleArray Phase(Nx,Ny,Nz);
|
||||
|
||||
// Extra copies of phi needed to compute time derivatives on CPU
|
||||
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;
|
||||
cDen = new double[2*N];
|
||||
cDistEven = new double[10*N];
|
||||
cDistOdd = new double[9*N];
|
||||
|
||||
IntArray LocalBlobID(Nx,Ny,Nz);
|
||||
DoubleArray Press(Nx,Ny,Nz);
|
||||
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);
|
||||
DoubleArray Vel_x(Nx,Ny,Nz); // Gradient of the phase indicator field
|
||||
DoubleArray Vel_y(Nx,Ny,Nz);
|
||||
DoubleArray Vel_z(Nx,Ny,Nz);
|
||||
/*****************************************************************
|
||||
VARIABLES FOR THE PMMC ALGORITHM
|
||||
****************************************************************** */
|
||||
//...........................................................................
|
||||
// Averaging variables
|
||||
//...........................................................................
|
||||
// local averages (to each MPI process)
|
||||
double trimdist=1.0; // pixel distance to trim surface for specified averages
|
||||
double awn,ans,aws,lwns,nwp_volume;
|
||||
double As, dummy;
|
||||
double vol_w, vol_n; // volumes the exclude the interfacial region
|
||||
double sat_w, sat_w_previous;
|
||||
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
|
||||
double Kwn,Kwn_global; // average Gaussian curavture - wn interface
|
||||
double KNwns,KNwns_global; // wns common curve normal curavture
|
||||
double KGwns,KGwns_global; // wns common curve geodesic curavture
|
||||
double trawn,trawn_global; // trimmed interfacial area
|
||||
double trJwn,trJwn_global; // trimmed interfacial area
|
||||
double trRwn,trRwn_global; // trimmed interfacial area
|
||||
DoubleArray van(3);
|
||||
DoubleArray vaw(3);
|
||||
DoubleArray vawn(3);
|
||||
DoubleArray vawns(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 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 vawns_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
|
||||
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 Values(20);
|
||||
DoubleArray KGwns_values(20);
|
||||
DoubleArray KNwns_values(20);
|
||||
DoubleArray ContactAngle(20);
|
||||
DoubleArray Curvature(20);
|
||||
DoubleArray DistValues(20);
|
||||
DoubleArray InterfaceSpeed(20);
|
||||
DoubleArray NormalVector(60);
|
||||
|
||||
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;
|
||||
|
||||
Point A,B,C,P;
|
||||
|
||||
// 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 n_nw_tris_beg, n_ns_tris_beg, n_ws_tris_beg;
|
||||
int c;
|
||||
//int newton_steps = 0;
|
||||
// Copy signed distance for device initialization
|
||||
CopyToDevice(dvcSignDist, Averages.SDs.data, dist_mem_size);
|
||||
//...........................................................................
|
||||
int ncubes = (Nx-2)*(Ny-2)*(Nz-2); // Exclude the "upper" halo
|
||||
IntArray cubeList(3,ncubes);
|
||||
int nc=0;
|
||||
//...........................................................................
|
||||
// Set up the cube list (very regular in this case due to lack of blob-ID)
|
||||
for (k=kstart; k<kfinish; k++){
|
||||
for (j=1; j<Ny-1; j++){
|
||||
for (i=1; i<Nx-1; i++){
|
||||
cubeList(0,nc) = i;
|
||||
cubeList(1,nc) = j;
|
||||
cubeList(2,nc) = k;
|
||||
nc++;
|
||||
}
|
||||
}
|
||||
}
|
||||
ncubes = nc;
|
||||
|
||||
int logcount = 0; // number of surface write-outs
|
||||
|
||||
@ -1127,84 +1002,15 @@ int main(int argc, char **argv)
|
||||
//......................................................................
|
||||
InitD3Q7(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
|
||||
InitD3Q7(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
|
||||
//.......................................................................
|
||||
// Once phase has been initialized, map solid to account for 'smeared' interface
|
||||
//......................................................................
|
||||
for (i=0; i<N; i++) SignDist.data[i] -= (1.0); //
|
||||
for (i=0; i<N; i++) Averages.SDs.data[i] -= (1.0); //
|
||||
//.......................................................................
|
||||
// Finalize setup for averaging domain
|
||||
Averages.SetupCubes(Dm);
|
||||
Averages.UpdateSolid();
|
||||
//.......................................................................
|
||||
|
||||
//...........................................................................
|
||||
// Gradient of the Signed Distance function
|
||||
//...........................................................................
|
||||
pmmc_MeshGradient(SignDist,SignDist_x,SignDist_y,SignDist_z,Nx,Ny,Nz);
|
||||
//...........................................................................
|
||||
CommunicateMeshHalo(SignDist_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(SignDist_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(SignDist_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 the phase indicator field and reset Copy, Den
|
||||
//*************************************************************************
|
||||
@ -1318,7 +1124,7 @@ int main(int argc, char **argv)
|
||||
//...........................................................................
|
||||
// Copy the phase indicator field for the earlier timestep
|
||||
DeviceBarrier();
|
||||
CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
|
||||
CopyToHost(Averages.Phase_tplus.data,Phi,N*sizeof(double));
|
||||
//...........................................................................
|
||||
//...........................................................................
|
||||
// Copy the data for for the analysis timestep
|
||||
@ -1327,14 +1133,15 @@ int main(int argc, char **argv)
|
||||
//...........................................................................
|
||||
DeviceBarrier();
|
||||
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
CopyToHost(Phase.data,Phi,N*sizeof(double));
|
||||
CopyToHost(Press.data,Pressure,N*sizeof(double));
|
||||
CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double));
|
||||
CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double));
|
||||
CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double));
|
||||
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
|
||||
CopyToHost(Averages.Press.data,Pressure,N*sizeof(double));
|
||||
CopyToHost(Averages.Vel_x.data,&Velocity[0],N*sizeof(double));
|
||||
CopyToHost(Averages.Vel_y.data,&Velocity[N],N*sizeof(double));
|
||||
CopyToHost(Averages.Vel_z.data,&Velocity[2*N],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);
|
||||
@ -1346,26 +1153,9 @@ int main(int argc, char **argv)
|
||||
//.........................................
|
||||
|
||||
sendtag = recvtag = 5;
|
||||
FILE *TIMELOG;
|
||||
if (rank==0){
|
||||
TIMELOG= fopen("timelog.tcat","a+");
|
||||
if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR)){
|
||||
// If timelog is empty, write a short header to list the averages
|
||||
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
|
||||
fprintf(TIMELOG,"time dEs "); // Timestep, Change in Surface Energy
|
||||
fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns sgkvpmawns KNwns KGwns "); // Scalar averages
|
||||
fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages
|
||||
fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz ");
|
||||
fprintf(TIMELOG,"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
|
||||
fprintf(TIMELOG,"Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
|
||||
fprintf(TIMELOG,"Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
|
||||
fprintf(TIMELOG,"trawn trJwn trRwn\n"); // trimmed curvature for wn surface
|
||||
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
|
||||
}
|
||||
}
|
||||
|
||||
err = 1.0;
|
||||
sat_w_previous = 1.01; // slightly impossible value!
|
||||
double sat_w_previous = 1.01; // slightly impossible value!
|
||||
if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
|
||||
//************ MAIN ITERATION LOOP ***************************************/
|
||||
while (timestep < timestepMax && err > tol ){
|
||||
@ -1758,7 +1548,8 @@ int main(int argc, char **argv)
|
||||
//...........................................................................
|
||||
// Copy the phase indicator field for the earlier timestep
|
||||
DeviceBarrier();
|
||||
CopyToHost(Phase_tplus.data,Phi,N*sizeof(double));
|
||||
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
|
||||
Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.Phase_tplus.data);
|
||||
//...........................................................................
|
||||
}
|
||||
if (timestep%1000 == 0){
|
||||
@ -1769,494 +1560,29 @@ int main(int argc, char **argv)
|
||||
//...........................................................................
|
||||
DeviceBarrier();
|
||||
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
CopyToHost(Phase.data,Phi,N*sizeof(double));
|
||||
CopyToHost(Press.data,Pressure,N*sizeof(double));
|
||||
CopyToHost(Vel_x.data,&Velocity[0],N*sizeof(double));
|
||||
CopyToHost(Vel_y.data,&Velocity[N],N*sizeof(double));
|
||||
CopyToHost(Vel_z.data,&Velocity[2*N],N*sizeof(double));
|
||||
CopyToHost(Averages.Phase.data,Phi,N*sizeof(double));
|
||||
CopyToHost(Averages.Press.data,Pressure,N*sizeof(double));
|
||||
CopyToHost(Averages.Vel_x.data,&Velocity[0],N*sizeof(double));
|
||||
CopyToHost(Averages.Vel_y.data,&Velocity[N],N*sizeof(double));
|
||||
CopyToHost(Averages.Vel_z.data,&Velocity[2*N],N*sizeof(double));
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
}
|
||||
if (timestep%1000 == 5){
|
||||
//...........................................................................
|
||||
// Copy the phase indicator field for the later timestep
|
||||
DeviceBarrier();
|
||||
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));
|
||||
//...........................................................................
|
||||
|
||||
//...........................................................................
|
||||
CommunicateMeshHalo(Phase, 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 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
|
||||
//...........................................................................
|
||||
// Pressure
|
||||
//...........................................................................
|
||||
CommunicateMeshHalo(Vel_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);
|
||||
//...........................................................................
|
||||
// Velocity
|
||||
//...........................................................................
|
||||
CommunicateMeshHalo(Press, 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(Vel_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(Vel_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);
|
||||
//...........................................................................
|
||||
// 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)
|
||||
//...........................................................................
|
||||
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;
|
||||
van(0) = van(1) = van(2) = 0.0;
|
||||
vawn(0) = vawn(1) = vawn(2) = 0.0;
|
||||
vawns(0) = vawns(1) = vawns(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;
|
||||
KGwns = KNwns = 0.0;
|
||||
Jwn = Kwn = efawns = 0.0;
|
||||
trJwn = trawn = trRwn = 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);
|
||||
|
||||
//...........................................................................
|
||||
// Compute volume averages
|
||||
for (int p=0;p<8;p++){
|
||||
|
||||
double delphi;
|
||||
if ( SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
|
||||
// 1-D index for this cube corner
|
||||
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
|
||||
// compute the norm of the gradient of the phase indicator field
|
||||
delphi = sqrt(Phase_x.data[n]*Phase_x.data[n]+Phase_y.data[n]*Phase_y.data[n]+Phase_z.data[n]*Phase_z.data[n]);
|
||||
// Compute the non-wetting phase volume contribution
|
||||
if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
|
||||
nwp_volume += 0.125;
|
||||
// volume the excludes the interfacial region
|
||||
if (delphi < 1e-4){
|
||||
vol_n += 0.125;
|
||||
// pressure
|
||||
pan += 0.125*Press.data[n];
|
||||
// velocity
|
||||
van(0) += 0.125*Vel_x.data[n];
|
||||
van(1) += 0.125*Vel_y.data[n];
|
||||
van(2) += 0.125*Vel_z.data[n];
|
||||
}
|
||||
}
|
||||
else if (delphi < 1e-4){
|
||||
// volume the excludes the interfacial region
|
||||
vol_w += 0.125;
|
||||
// pressure
|
||||
paw += 0.125*Press.data[n];
|
||||
// velocity
|
||||
vaw(0) += 0.125*Vel_x.data[n];
|
||||
vaw(1) += 0.125*Vel_y.data[n];
|
||||
vaw(2) += 0.125*Vel_z.data[n];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//...........................................................................
|
||||
// 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);
|
||||
|
||||
// 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);
|
||||
|
||||
// Integrate the mean curvature
|
||||
Jwn += pmmc_CubeSurfaceInterpValue(CubeValues,MeanCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
|
||||
Kwn += pmmc_CubeSurfaceInterpValue(CubeValues,GaussCurvature,nw_pts,nw_tris,Values,i,j,k,n_nw_pts,n_nw_tris);
|
||||
|
||||
// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
|
||||
pmmc_CubeTrimSurfaceInterpValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues,
|
||||
i,j,k,n_nw_pts,n_nw_tris,trimdist,trawn,trJwn);
|
||||
|
||||
pmmc_CubeTrimSurfaceInterpInverseValues(CubeValues,MeanCurvature,SignDist,nw_pts,nw_tris,Values,DistValues,
|
||||
i,j,k,n_nw_pts,n_nw_tris,trimdist,dummy,trRwn);
|
||||
|
||||
// Compute the normal speed of the 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_CommonCurveSpeed(CubeValues, dPdt, vawns,Phase_x,Phase_y,Phase_z,SignDist_x,SignDist_y,SignDist_z,
|
||||
local_nws_pts,i,j,k,n_local_nws_pts);
|
||||
|
||||
pmmc_CurveCurvature(Phase, SignDist, KNwns_values, KGwns_values, KNwns, KGwns,
|
||||
nws_pts, n_nws_pts, i, j, k);
|
||||
|
||||
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);
|
||||
MPI_Allreduce(&nwp_volume,&nwp_volume_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
MPI_Allreduce(&awn,&awn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
MPI_Allreduce(&ans,&ans_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
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(&Kwn,&Kwn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
MPI_Allreduce(&KGwns,&KGwns_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
MPI_Allreduce(&KNwns,&KNwns_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(&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(&vawns(0),&vawns_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_Allreduce(&trawn,&trawn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
MPI_Allreduce(&trJwn,&trJwn_global,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
|
||||
MPI_Allreduce(&trRwn,&trRwn_global,1,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
|
||||
// See McClure, Prins and Miller (2014)
|
||||
//.........................................................................
|
||||
dAwn += awn_global;
|
||||
dAns += ans_global;
|
||||
dEs = 6.01603*alpha*(dAwn + dAns);
|
||||
dAwn = -awn_global; // Get ready for the next analysis interval
|
||||
dAns = -ans_global;
|
||||
|
||||
// Normalize the phase averages
|
||||
// (density of both components = 1.0)
|
||||
if (vol_w_global > 0.0){
|
||||
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;
|
||||
}
|
||||
if (vol_n_global > 0.0){
|
||||
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
|
||||
if (awn_global > 0.0){
|
||||
Jwn_global /= awn_global;
|
||||
Kwn_global /= awn_global;
|
||||
for (i=0; i<3; i++) vawn_global(i) /= awn_global;
|
||||
for (i=0; i<6; i++) Gwn_global(i) /= awn_global;
|
||||
}
|
||||
if (lwns_global > 0.0){
|
||||
efawns_global /= lwns_global;
|
||||
KNwns_global /= lwns_global;
|
||||
KGwns_global /= lwns_global;
|
||||
for (i=0; i<3; i++) vawns_global(i) /= lwns_global;
|
||||
}
|
||||
if (trawn_global > 0.0){
|
||||
trJwn_global /= trawn_global;
|
||||
trRwn_global /= trawn_global;
|
||||
trRwn_global = 2.0*fabs(trRwn_global);
|
||||
trJwn_global = fabs(trJwn_global);
|
||||
}
|
||||
|
||||
if (ans_global > 0.0) for (i=0; i<6; i++) Gns_global(i) /= ans_global;
|
||||
if (aws_global > 0.0) for (i=0; i<6; i++) Gws_global(i) /= aws_global;
|
||||
|
||||
D = 6.0*(1-porosity)/(ans_global+aws_global)/iVol_global;
|
||||
MPI_Bcast(&D,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
|
||||
MPI_Barrier(MPI_COMM_WORLD);
|
||||
|
||||
//sat_w = 1.0 - nwp_volume_global*iVol_global/porosity;
|
||||
sat_w = 1.0 - nwp_volume_global/pore_vol;
|
||||
// Compute the specific interfacial areas and common line length (dimensionless per unit volume)
|
||||
awn_global = awn_global*iVol_global*D;
|
||||
ans_global = ans_global*iVol_global*D;
|
||||
aws_global = aws_global*iVol_global*D;
|
||||
dEs = dEs*iVol_global*D;
|
||||
lwns_global = lwns_global*iVol_global*D*D;
|
||||
|
||||
//.........................................................................
|
||||
if (rank==0){
|
||||
printf("Sauter Mean Diameter = %f \n",D);
|
||||
fprintf(TIMELOG,"%i %.5g ",timestep-5,dEs); // change in surface energy
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas
|
||||
fprintf(TIMELOG,"%.5g %.5g ",Jwn_global, Kwn_global); // curvature of wn interface
|
||||
fprintf(TIMELOG,"%.5g ",lwns_global); // common curve length
|
||||
fprintf(TIMELOG,"%.5g ",efawns_global); // average contact angle
|
||||
fprintf(TIMELOG,"%.5g %.5g ",KNwns_global, KGwns_global); // curvature of wn interface
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g ",vawns_global(0),vawns_global(1),vawns_global(2)); // velocity of wn interface
|
||||
fprintf(TIMELOG,"%.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
|
||||
fprintf(TIMELOG,"%.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
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",
|
||||
Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface
|
||||
fprintf(TIMELOG,"%.5g %.5g %.5g\n",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature
|
||||
fflush(TIMELOG);
|
||||
}
|
||||
CopyToHost(Averages.Phase_tminus.data,Phi,N*sizeof(double));
|
||||
Averages.ColorToSignedDistance(beta,Averages.Phase_tminus.data,Averages.Phase_tminus.data);
|
||||
//....................................................................
|
||||
// The following need to be called each time new averages are computed
|
||||
Averages.Initialize();
|
||||
Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.SDn.data);
|
||||
Averages.UpdateMeshValues();
|
||||
Averages.ComputeLocal();
|
||||
Averages.Reduce();
|
||||
Averages.PrintAll(timestep);
|
||||
//....................................................................
|
||||
}
|
||||
|
||||
if (timestep%RESTART_INTERVAL == 0){
|
||||
if (pBC){
|
||||
//err = fabs(sat_w - sat_w_previous);
|
||||
@ -2412,68 +1738,20 @@ int main(int argc, char **argv)
|
||||
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
|
||||
if (rank==0) printf("********************************************************\n");
|
||||
|
||||
// printf("Local File Name = %s \n",LocalRankFilename);
|
||||
FILE *FINALSTATE;
|
||||
if (rank==0){
|
||||
fclose(TIMELOG);
|
||||
FINALSTATE= fopen("finalstate.tcat","a");
|
||||
fprintf(FINALSTATE,"%i %.5g ",timestep-5,dEs); // change in surface energy
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas
|
||||
fprintf(FINALSTATE,"%.5g %5g ",Jwn_global, Kwn_global); // curvature of wn interface
|
||||
fprintf(FINALSTATE,"%.5g ",lwns_global); // common curve length
|
||||
fprintf(FINALSTATE,"%.5g ",efawns_global); // average contact angle
|
||||
fprintf(FINALSTATE,"%.5g %5g ",KNwns_global, KGwns_global); // curvature of wns curve
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g ",vaw_global(0),vaw_global(1),vaw_global(2)); // average velocity of w phase
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g ",van_global(0),van_global(1),van_global(2)); // average velocity of n phase
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g ",vawn_global(0),vawn_global(1),vawn_global(2)); // velocity of wn interface
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g ",vawns_global(0),vawns_global(1),vawns_global(2)); // velocity of wn interface
|
||||
fprintf(FINALSTATE,"%.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
|
||||
fprintf(FINALSTATE,"%.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
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g %.5g %.5g %.5g ",
|
||||
Gws_global(0),Gws_global(1),Gws_global(2),Gws_global(3),Gws_global(4),Gws_global(5)); // orientation of ws interface
|
||||
fprintf(FINALSTATE,"%.5g %.5g %.5g\n",trawn_global, trJwn_global, trRwn_global); // Trimmed curvature
|
||||
fclose(FINALSTATE);
|
||||
}
|
||||
|
||||
|
||||
//#ifdef WriteOutput
|
||||
CopyToHost(Phase.data,Phi,N*sizeof(double));
|
||||
sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
|
||||
FILE *PHASE;
|
||||
PHASE = fopen(LocalRankFilename,"wb");
|
||||
fwrite(Phase.data,8,N,PHASE);
|
||||
fwrite(Averages.Phase.data,8,N,PHASE);
|
||||
// fwrite(MeanCurvature.data,8,N,PHASE);
|
||||
fclose(PHASE);
|
||||
//#endif
|
||||
ComputePressureD3Q19(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
CopyToHost(Press.data,Pressure,N*sizeof(double));
|
||||
CommunicateMeshHalo(Press, 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);
|
||||
|
||||
sprintf(LocalRankFilename,"%s%s","Pressure.",LocalRankString);
|
||||
FILE *PRESS;
|
||||
PRESS = fopen(LocalRankFilename,"wb");
|
||||
fwrite(Press.data,8,N,PRESS);
|
||||
fwrite(Averages.Press.data,8,N,PRESS);
|
||||
fclose(PRESS);
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user