Merge branch 'morphLBM' of github.com:JamesEMcClure/LBPM-WIA into morphLBM
This commit is contained in:
JamesEMcclure
commit
a6ade4eb0e
|
|
@ -564,7 +564,6 @@ void ScaLBL_ColorModel::Run(){
|
|||
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_DeviceBarrier();
|
||||
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
|
||||
|
|
@ -649,7 +648,6 @@ void ScaLBL_ColorModel::Run(){
|
|||
double volA = Averages->gnb.V;
|
||||
volA /= Dm->Volume;
|
||||
volB /= Dm->Volume;;
|
||||
initial_volume = volA*Dm->Volume;
|
||||
double vA_x = Averages->gnb.Px/Averages->gnb.M;
|
||||
double vA_y = Averages->gnb.Py/Averages->gnb.M;
|
||||
double vA_z = Averages->gnb.Pz/Averages->gnb.M;
|
||||
|
|
@ -683,9 +681,10 @@ void ScaLBL_ColorModel::Run(){
|
|||
isSteady = true;
|
||||
|
||||
if ( isSteady ){
|
||||
initial_volume = volA*Dm->Volume;
|
||||
MORPH_ADAPT = true;
|
||||
CURRENT_MORPH_TIMESTEPS=0;
|
||||
delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
|
||||
delta_volume_target = Dm->Volume*morph_delta; //*volA ???? // set target volume change
|
||||
Averages->Full();
|
||||
Averages->Write(timestep);
|
||||
analysis.WriteVisData( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
|
||||
|
|
@ -764,7 +763,8 @@ void ScaLBL_ColorModel::Run(){
|
|||
else if (USE_SEED){
|
||||
delta_volume = volA*Dm->Volume - initial_volume;
|
||||
CURRENT_MORPH_TIMESTEPS += analysis_interval;
|
||||
double massChange = SeedPhaseField(seed_water);
|
||||
double effective_seed_water = seed_water*(1.0+volB/volA);
|
||||
double massChange = SeedPhaseField(effective_seed_water);
|
||||
if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target);
|
||||
}
|
||||
else if (USE_MORPHOPEN_OIL){
|
||||
|
|
@ -983,10 +983,16 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
|
|||
srand(time(NULL));
|
||||
double mass_loss =0.f;
|
||||
double count =0.f;
|
||||
DoubleArray phase(Nx,Ny,Nz);
|
||||
double *Aq_tmp, *Bq_tmp;
|
||||
double SEED_THRESHOLD = -0.9;
|
||||
|
||||
Aq_tmp = new double [7*Np];
|
||||
Bq_tmp = new double [7*Np];
|
||||
|
||||
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double));
|
||||
ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double));
|
||||
|
||||
/* for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
double random_value = double(rand())/ RAND_MAX;
|
||||
|
|
@ -1005,26 +1011,70 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
|
|||
}
|
||||
}
|
||||
}
|
||||
*/
|
||||
double oil_value = 0.0;
|
||||
double water_value = 1.0;
|
||||
for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
|
||||
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
|
||||
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
|
||||
double phase_id = (dA - dB) / (dA + dB);
|
||||
double random_value = double(rand())/ RAND_MAX;
|
||||
if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
|
||||
Aq_tmp[n] = 0.3333333333333333*oil_value;
|
||||
Aq_tmp[n+Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value;
|
||||
|
||||
Bq_tmp[n] = 0.3333333333333333*water_value;
|
||||
Bq_tmp[n+Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+2*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+3*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+4*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+5*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+6*Np] = 0.1111111111111111*water_value;
|
||||
mass_loss += oil_value - dA;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){
|
||||
double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
|
||||
double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
|
||||
double phase_id = (dA - dB) / (dA + dB);
|
||||
double random_value = double(rand())/ RAND_MAX;
|
||||
if (phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){
|
||||
Aq_tmp[n] = 0.3333333333333333*oil_value;
|
||||
Aq_tmp[n+Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value;
|
||||
Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value;
|
||||
|
||||
Bq_tmp[n] = 0.3333333333333333*water_value;
|
||||
Bq_tmp[n+Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+2*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+3*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+4*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+5*Np] = 0.1111111111111111*water_value;
|
||||
Bq_tmp[n+6*Np] = 0.1111111111111111*water_value;
|
||||
mass_loss += oil_value - dA;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
count= sumReduce( Dm->Comm, count);
|
||||
mass_loss= sumReduce( Dm->Comm, mass_loss);
|
||||
if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
|
||||
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
|
||||
|
||||
// 7. Re-initialize phase field and density
|
||||
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){
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
|
||||
}
|
||||
if (Dm->kproc() == nprocz-1){
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
|
||||
}
|
||||
}
|
||||
// Need to initialize Aq, Bq, Den, Phi directly
|
||||
//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
|
||||
ScaLBL_CopyToDevice(Aq, Aq_tmp, 7*Np*sizeof(double));
|
||||
ScaLBL_CopyToDevice(Bq, Bq_tmp, 7*Np*sizeof(double));
|
||||
|
||||
return(mass_loss);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -38,7 +38,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;
|
||||
|
||||
|
|
@ -218,12 +218,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
|
||||
|
|
|
|||
|
|
@ -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
234
tests/TestTopo3D.cpp
Normal 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;
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user