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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure 2018-07-28 16:03:51 -04:00
commit 4063162b7d
6 changed files with 83 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
analysis/Minkowski.cpp
View File

@ -67,6 +67,22 @@ void Minkowski::Initialize()
vol_n = euler = Jn = An = Kn = 0.0;
}
double Minkowski::V(){
return vol_n_global;
}
double Minkowski::A(){
return An_global;
}
double Minkowski::J(){
return Jn_global;
}
double Minkowski::X(){
return euler_global_global;
}
void Minkowski::UpdateMeshValues()
{

6
analysis/Minkowski.h
View File

@ -69,6 +69,11 @@ public:
DoubleArray SDn_y;
DoubleArray SDn_z;
double V();
double A();
double J();
double X();
//...........................................................................
Minkowski(std::shared_ptr <Domain> Dm);
~Minkowski();
@ -80,6 +85,7 @@ public:
void PrintAll();
int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel);
void SortBlobs();
};

2
common/ScaLBL.cpp
View File

@ -1121,14 +1121,12 @@ void ScaLBL_Communicator::RecvGrad(double *phi, double *grad){
ScaLBL_DeviceBarrier();
//...................................................................................
// NOTE: AA Routine writes to opposite
// Unpack the gradributions on the device
//...................................................................................
//...Unpacking for x face(2,8,10,12,14)................................
ScaLBL_Gradient_Unpack(1.0,-1,0,0,dvcRecvDist_x,0,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,-1,0,dvcRecvDist_x,recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,1,0,dvcRecvDist_x,2*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,0,-1,dvcRecvDist_x,3*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
ScaLBL_Gradient_Unpack(0.5,-1,0,1,dvcRecvDist_x,4*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
//...................................................................................
//...Packing for X face(1,7,9,11,13)................................

24
models/ColorModel.cpp
View File

@ -246,8 +246,8 @@ void ScaLBL_ColorModel::Create(){
}
}
// check that TmpMap is valid
for (int idx=0; idx<ScaLBL_Comm->last_interior; idx++){
if (idx == ScaLBL_Comm->next) idx = ScaLBL_Comm->first_interior;
for (int idx=0; idx<ScaLBL_Comm->LastInterior(); idx++){
if (idx == ScaLBL_Comm->LastExterior()) idx = ScaLBL_Comm->FirstInterior();
int n = TmpMap[idx];
if (n > Nx*Ny*Nz){
printf("Bad value! idx=%i \n");
@ -316,8 +316,8 @@ void ScaLBL_ColorModel::Initialize(){
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf ("Initializing phase field \n");
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->next, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (BoundaryCondition >0 ){
if (Dm->kproc()==0){
@ -365,9 +365,9 @@ void ScaLBL_ColorModel::Run(){
// Compute the Phase indicator field
// Read for Aq, Bq happens in this routine (requires communication)
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
@ -375,7 +375,7 @@ void ScaLBL_ColorModel::Run(){
ScaLBL_Comm_Regular->SendHalo(Phi);
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set BCs
@ -392,23 +392,23 @@ void ScaLBL_ColorModel::Run(){
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->next, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
// Set boundary conditions
@ -425,7 +425,7 @@ void ScaLBL_ColorModel::Run(){
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->next, Np);
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************

55
models/MRTModel.cpp
View File

@ -124,6 +124,8 @@ void ScaLBL_MRTModel::Create(){
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
Morphology = new Minkowski(Mask);
}
void ScaLBL_MRTModel::Initialize(){
@ -162,21 +164,56 @@ void ScaLBL_MRTModel::Run(){
}
//************************************************************************/
stoptime = MPI_Wtime();
// cout << "CPU time: " << (stoptime - starttime) << " seconds" << endl;
cputime = stoptime - starttime;
// cout << "Lattice update rate: "<< double(Nx*Ny*Nz*timestep)/cputime/1000000 << " MLUPS" << endl;
double MLUPS = double(Np*timestep)/cputime/1000000;
// if (rank==0) printf("********************************************************\n");
// if (rank==0) printf("CPU time = %f \n", cputime);
// if (rank==0) printf("Lattice update rate (per process)= %f MLUPS \n", MLUPS);
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
// Performance obtained from each node
double MLUPS = double(Np)/cputime/1000000;
if (rank==0) printf("********************************************************\n");
if (rank==0) printf("CPU time = %f \n", cputime);
if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
MLUPS *= nprocs;
if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
if (rank==0) printf("********************************************************\n");
}
void ScaLBL_MRTModel::VelocityField(double *Vz){
void ScaLBL_MRTModel::VelocityField(double *VELOCITY){
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_CopyToHost(&Vz[0],&Velocity[2*Np],SIZE);
ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],SIZE);
memcpy(Morphology.SDn.data(), Distance.data(), N*sizeof(double));
Morphology.Initialize();
Morphology.UpdateMeshValues();
Morphology.ComputeLocal();
Morphology.Reduce();
unsigned long int count_loc=0;
unsigned long int count;
double vax,vay,vaz;
double vax_loc,vay_loc,vaz_loc;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int n=ScaLBL_Comm->0; n<ScaLBL_Comm->LastExterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
}
for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
}
MPI_Allreduce(vax_loc,vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(vay_loc,vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(vaz_loc,vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
if (rank==0) printf("Vs As Js Xs vx vy vz\n");
if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
}

3
models/MRTModel.h
View File

@ -12,6 +12,7 @@
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
class ScaLBL_MRTModel{
@ -54,6 +55,8 @@ public:
double *fq;
double *Velocity;
double *Pressure;
Minkowski Morphology;
private:
MPI_Comm comm;