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

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

Browse Source
This commit is contained in:
James E McClure 2019-06-11 15:28:38 -04:00
commit b565fa7820
6 changed files with 245 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
analysis/SubPhase.cpp
View File

@ -245,8 +245,8 @@ void SubPhase::Basic(){
force_mag = 1.0;
}
double saturation=gwb.V/(gwb.V + gnb.V);
double water_flow_rate=gwb.V*sqrt(gwb.Px*gwb.Px + gwb.Py*gwb.Py + gwb.Pz*gwb.Pz)/gwb.M / Dm->Volume;
double not_water_flow_rate=gnb.V*sqrt(gnb.Px*gnb.Px + gnb.Py*gnb.Py + gnb.Pz*gnb.Pz)/gnb.M/ Dm->Volume;
double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume;
double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume;
double total_flow_rate = water_flow_rate + not_water_flow_rate;
double fractional_flow= water_flow_rate / total_flow_rate;

8
gpu/D3Q19.cu
View File

@ -2499,9 +2499,8 @@ extern "C" double ScaLBL_D3Q19_AAeven_Flux_BC_z(int *list, double *dist, double
// Allocate memory to store the sums
double din;
double *sum;
double sum[1];
double *dvcsum;
cudaMallocHost((void **)&sum,sizeof(double));
cudaMalloc((void **)&dvcsum,sizeof(double)*count);
cudaMemset(dvcsum,0,sizeof(double)*count);
int sharedBytes = 512*sizeof(double);
@ -2544,9 +2543,8 @@ extern "C" double ScaLBL_D3Q19_AAodd_Flux_BC_z(int *neighborList, int *list, dou
// Allocate memory to store the sums
double din;
double *sum;
double sum[1];
double *dvcsum;
cudaMallocHost((void **)&sum,sizeof(double));
cudaMalloc((void **)&dvcsum,sizeof(double)*count);
cudaMemset(dvcsum,0,sizeof(double)*count);
int sharedBytes = 512*sizeof(double);
@ -2595,7 +2593,7 @@ extern "C" double ScaLBL_D3Q19_Flux_BC_Z(double *disteven, double *distodd, doub
dvc_D3Q19_Flux_BC_Z<<<GRID,512>>>(disteven, distodd, flux, dvcsum, Nx, Ny, Nz, outlet);
// Now read the total flux
cudaMemcpy(&sum[0],&dvcsum[0],sizeof(double),cudaMemcpyDeviceToHost);
cudaMemcpy(&sum[0],dvcsum,sizeof(double),cudaMemcpyDeviceToHost);
// free the memory needed for reduction

4
models/ColorModel.cpp
View File

@ -654,8 +654,8 @@ void ScaLBL_ColorModel::Run(){
force_mag = 1.0;
}
double current_saturation = volB/(volA+volB);
double flow_rate_A = volA*sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
double flow_rate_B = volB*sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
if ( morph_timesteps > morph_interval ){

6
models/MRTModel.cpp
View File

@ -23,7 +23,7 @@ void ScaLBL_MRTModel::ReadParams(string filename){
tau = 1.0;
timestepMax = 100000;
tolerance = 0.01;
tolerance = 1.0e-8;
Fx = Fy = 0.0;
Fz = 1.0e-5;
@ -203,12 +203,12 @@ void ScaLBL_MRTModel::Run(){
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
starttime = MPI_Wtime();
if (rank==0) printf("Beginning AA timesteps...\n");
if (rank==0) printf("Beginning AA timesteps, timestepMax = %i \n", timestepMax);
if (rank==0) printf("********************************************************\n");
timestep=0;
double error = 1.0;
double flow_rate_previous = 0.0;
while (timestep < timestepMax && error < tolerance) {
while (timestep < timestepMax && error > tolerance) {
//************************************************************************/
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL

1
tests/CMakeLists.txt
View File

@ -41,6 +41,7 @@ CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_
ADD_LBPM_TEST( pmmc_cylinder )
ADD_LBPM_TEST( TestTorus )
ADD_LBPM_TEST( TestTorusEvolve )
ADD_LBPM_TEST( TestTopo3D )
ADD_LBPM_TEST( TestFluxBC )
ADD_LBPM_TEST( TestMap )
#ADD_LBPM_TEST( TestMRT )

234
tests/TestTopo3D.cpp Normal file
View File

@ -0,0 +1,234 @@
// Sequential blob analysis
// Reads parallel simulation data and performs connectivity analysis
// and averaging on a blob-by-blob basis
// James E. McClure 2014
#include <iostream>
#include <math.h>
#include "common/Communication.h"
#include "analysis/analysis.h"
#include "analysis/Minkowski.h"
#include "IO/MeshDatabase.h"
std::shared_ptr<Database> loadInputs( int nprocs )
{
//auto db = std::make_shared<Database>( "Domain.in" );
auto db = std::make_shared<Database>();
db->putScalar<int>( "BC", 0 );
db->putVector<int>( "nproc", { 1, 1, 1 } );
db->putVector<int>( "n", { 100, 100, 100 } );
db->putScalar<int>( "nspheres", 1 );
db->putVector<double>( "L", { 1, 1, 1 } );
return db;
}
int main(int argc, char **argv)
{
// Initialize MPI
int rank, nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
{ // Limit scope so variables that contain communicators will free before MPI_Finialize
if ( rank==0 ) {
printf("-----------------------------------------------------------\n");
printf("Unit test 3D topologies \n");
printf("-----------------------------------------------------------\n");
}
//.......................................................................
// Reading the domain information file
//.......................................................................
int i,j,k,n;
// Load inputs
auto db = loadInputs( nprocs );
int Nx = db->getVector<int>( "n" )[0];
int Ny = db->getVector<int>( "n" )[1];
int Nz = db->getVector<int>( "n" )[2];
int nprocx = db->getVector<int>( "nproc" )[0];
int nprocy = db->getVector<int>( "nproc" )[1];
int nprocz = db->getVector<int>( "nproc" )[2];
if (rank==0){
printf("********************************************************\n");
printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
printf("********************************************************\n");
}
// Get the rank info
std::shared_ptr<Domain> Dm(new Domain(db,comm));
Nx += 2;
Ny += 2;
Nz += 2;
int N = Nx*Ny*Nz;
//.......................................................................
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
Dm->id[n] = 1;
}
}
}
//.......................................................................
Dm->CommInit(); // Initialize communications for domains
//.......................................................................
// Create visualization structure
std::vector<IO::MeshDataStruct> visData;
fillHalo<double> fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);;
IO::initialize("","silo","false");
// Create the MeshDataStruct
visData.resize(1);
visData[0].meshName = "domain";
visData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
auto PhaseVar = std::make_shared<IO::Variable>();
PhaseVar->name = "phase";
PhaseVar->type = IO::VariableType::VolumeVariable;
PhaseVar->dim = 1;
PhaseVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(PhaseVar);
//.......................................................................
// Assign the phase ID field based and the signed distance
//.......................................................................
double R1,R2,R;
double CX,CY,CZ; //CY1,CY2;
CX=Nx*nprocx*0.5;
CY=Ny*nprocy*0.5;
CZ=Nz*nprocz*0.5;
R1 = (Nx-2)*nprocx*0.3; // middle radius
R2 = (Nx-2)*nprocx*0.1; // donut thickness
R = 0.4*nprocx*(Nx-2);
Minkowski Object(Dm);
int timestep = 0;
double x,y,z;
// partial torus
timestep += 1;
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
// global position relative to center
x = Dm->iproc()*(Nx-2)+i - CX - 0.1;
y = Dm->jproc()*(Ny-2)+j - CY - 0.1;
z = Dm->kproc()*(Nz-2)+k - CZ -0.1;
//..............................................................................
if (x <= 0 || y<=0) {
// Single torus
Object.distance(i,j,k) = R2 - sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z);
}
else {
double d1 = R2-sqrt(x*x +(y-R1)*(y-R1) + z*z);
double d2 = R2-sqrt((x-R1)*(x-R1)+y*y + z*z);
Object.distance(i,j,k) = max(d1,d2);
}
if (Object.distance(i,j,k) > 0.0){
Dm->id[n] = 2;
Object.id(i,j,k) = 2;
}
else{
Dm->id[n] = 1;
Object.id(i,j,k) = 1;
}
}
}
}
ASSERT(visData[0].vars[0]->name=="phase");
Array<double>& PhaseData = visData[0].vars[0]->data;
fillData.copy(Object.distance,PhaseData);
IO::writeData( timestep, visData, comm );
//spherical shell
timestep += 1;
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
// global position relative to center
x = Dm->iproc()*(Nx-2)+i - CX - 0.1;
y = Dm->jproc()*(Ny-2)+j - CY - 0.1;
z = Dm->kproc()*(Nz-2)+k - CZ - 0.1;
//..............................................................................
// Single torus
double d1 = sqrt(x*x+y*y+z*z)-(R1-R2);
double d2 = R-sqrt(x*x+y*y+z*z);
Object.distance(i,j,k) = min(d1,d2);
if (Object.distance(i,j,k) > 0.0){
Dm->id[n] = 2;
Object.id(i,j,k) = 2;
}
else{
Dm->id[n] = 1;
Object.id(i,j,k) = 1;
}
}
}
}
ASSERT(visData[0].vars[0]->name=="phase");
PhaseData = visData[0].vars[0]->data;
fillData.copy(Object.distance,PhaseData);
IO::writeData( timestep, visData, comm );
// bowl
timestep += 1;
for ( k=1;k<Nz-1;k++){
for ( j=1;j<Ny-1;j++){
for ( i=1;i<Nx-1;i++){
n = k*Nx*Ny+j*Nx+i;
// global position relative to center
x = Dm->iproc()*(Nx-2)+i - CX - 0.1;
y = Dm->jproc()*(Ny-2)+j - CY - 0.1;
z = Dm->kproc()*(Nz-2)+k - CZ - 0.1;
//..............................................................................
// Bowl
if (z > 0 ){
Object.distance(i,j,k) = R2-sqrt((sqrt(x*x+y*y) - R1)*(sqrt(x*x+y*y) - R1) + z*z);
}
else
{
double d1 = sqrt(x*x+y*y+z*z)-(R1-R2);
double d2 = R-sqrt(x*x+y*y+z*z);
Object.distance(i,j,k) = min(d1,d2);
}
if (Object.distance(i,j,k) > 0.0){
Dm->id[n] = 2;
Object.id(i,j,k) = 2;
}
else{
Dm->id[n] = 1;
Object.id(i,j,k) = 1;
}
}
}
}
ASSERT(visData[0].vars[0]->name=="phase");
PhaseData = visData[0].vars[0]->data;
fillData.copy(Object.distance,PhaseData);
IO::writeData( timestep, visData, comm );
} // Limit scope so variables that contain communicators will free before MPI_Finialize
MPI_Barrier(comm);
MPI_Finalize();
return 0;
}