how doth the little crocodile

This commit is contained in:
James E McClure 2018-05-19 17:28:39 -04:00
parent b56cbcdf61
commit c18bf1792f
3 changed files with 31 additions and 42 deletions

View File

@ -40,6 +40,8 @@ void ScaLBL_ColorModel::ReadParams(string filename){
inletB=0.f;
outletA=0.f;
outletB=1.f;
if (BoundaryCondition==4) flux = din*rhoA; // mass flux must adjust for density (see formulation for details)
// Read domain parameters
auto L = domain_db->getVector<double>( "L" );
@ -56,19 +58,14 @@ void ScaLBL_ColorModel::ReadParams(string filename){
nprocy = nproc[1];
nprocz = nproc[2];
if (BoundaryCondition==4) flux = din*rhoA; // mass flux must adjust for density (see formulation for details)
}
void ScaLBL_ColorModel::SetDomain(){
Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis
Mask = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // mask domain removes immobile phases
Nx+=2; Ny+=2; Nz += 2;
N = Nx*Ny*Nz;
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
// local copy of the ids
id = new char[N];
MPI_Barrier(comm);
Dm->CommInit();
MPI_Barrier(comm);
@ -80,7 +77,7 @@ void ScaLBL_ColorModel::ReadInput(){
//.......................................................................
if (rank == 0) printf("Read input media... \n");
//.......................................................................
Dm->ReadIDs();
Mask->ReadIDs();
sprintf(LocalRankString,"%05d",Dm->rank());
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
@ -194,7 +191,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=Dm->id[n];
VALUE=Mask->id[n];
// Assign the affinity from the paired list
for (int idx=0; idx < NLABELS; idx++){
//printf("rank=%i, idx=%i, value=%i, %i, \n",rank(),idx, VALUE,LabelList[idx]);
@ -222,7 +219,7 @@ void ScaLBL_ColorModel::Create(){
//.........................................................
// Initialize communication structures in averaging domain
for (int i=0; i<Nx*Ny*Nz; i++) Mask->id[i] = Dm->id[i];
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
Mask->CommInit();
Np=Mask->PoreCount();
//...........................................................................
@ -285,13 +282,7 @@ void ScaLBL_ColorModel::Create(){
/********************************************************
* AssignComponentLabels *
********************************************************/
void ScaLBL_ColorModel::Initialize(){
/*
* This function initializes model (includes both mobile and immobile components)
*/
int rank=Dm->rank();
// Compute the solid interaction potential and copy result to device
void ScaLBL_ColorModel::SolidPotential(){
if (rank==0) printf("Computing solid interaction potential \n");
double *PhaseLabel;
PhaseLabel=new double [Nx*Ny*Nz];
@ -390,8 +381,15 @@ void ScaLBL_ColorModel::Initialize(){
PFILE = fopen(LocalRankFilename,"wb");
fwrite(Psnorm.data(),8,N,PFILE);
fclose(PFILE);
}
void ScaLBL_ColorModel::Initialize(){
/*
* This function initializes model
*/
int rank=Dm->rank();
double count_wet=0.f;
double *PhaseLabel;
PhaseLabel=new double [Nx*Ny*Nz];
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){

View File

@ -24,10 +24,11 @@ public:
// functions in they should be run
void ReadParams(string filename);
void ReadParams(std::shared_ptr<Database> db0);
void AssignComponentLabels(double *phase);
void SetDomain();
void ReadInput();
void Create();
void Initialize();
void SolidPotential();
void Run();
void WriteDebug();
@ -74,6 +75,7 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignComponentLabels(double *phase);
};

View File

@ -14,8 +14,6 @@ using namespace std;
inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
// initialize a bubble
int i,j,k,n;
int Np=0;
double sum=0.f;
int rank = ColorModel.Mask->rank();
int nprocx = ColorModel.Mask->nprocx();
int nprocy = ColorModel.Mask->nprocy();
@ -24,21 +22,11 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
int Ny = ColorModel.Mask->Ny;
int Nz = ColorModel.Mask->Nz;
if (rank == 0) cout << "Setting up bubble..." << endl;
sum=0; Np=0;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
ColorModel.Averages->SDs(i,j,k) = 100.f;
// Initialize phase positions field
if (ColorModel.Averages->SDs(i,j,k) < 0.0){
ColorModel.id[n] = 0;
ColorModel.id[n] = 0;
}
else {
sum++;
Np++;
}
}
}
}
@ -47,19 +35,19 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
int iglobal= i+(Nx-2)*ColorModel.Dm->iproc();
int jglobal= j+(Ny-2)*ColorModel.Dm->jproc();
int kglobal= k+(Nz-2)*ColorModel.Dm->kproc();
int iglobal= i+(Nx-2)*ColorModel.Mask->iproc();
int jglobal= j+(Ny-2)*ColorModel.Mask->jproc();
int kglobal= k+(Nz-2)*ColorModel.Mask->kproc();
// Initialize phase position field for parallel bubble test
if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
+(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
ColorModel.id[n] = 2;
ColorModel.id[n] = 2;
ColorModel.Mask->id[n] = 2;
ColorModel.Mask->id[n] = 2;
}
else{
ColorModel.id[n]=1;
ColorModel.id[n]=1;
ColorModel.Mask->id[n]=1;
ColorModel.Mask->id[n]=1;
}
}
}
@ -92,12 +80,13 @@ int main(int argc, char **argv)
auto filename = argv[1];
ScaLBL_ColorModel ColorModel(rank,nprocs,comm);
ColorModel.ReadParams(filename);
//ColorModel.ReadInput();
ColorModel.SetDomain();
//ColorModel.ReadInput();
double radius=15.5;
InitializeBubble(ColorModel,radius);
ColorModel.Create();
ColorModel.Initialize();
ColorModel.Run();
ColorModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
ColorModel.Run();
ColorModel.WriteDebug();
}
// ****************************************************