fix interface averaging for mass /momentum

This commit is contained in:
James E McClure 2019-03-26 10:26:02 -04:00
parent 2378a14553
commit ebabdcaef0
2 changed files with 33 additions and 16 deletions

View File

@ -290,7 +290,7 @@ void SubPhase::Full(){
if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y; if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z; if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z;
nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset(); nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset(); iwnc.reset();
Dm->CommunicateMeshHalo(Phi); Dm->CommunicateMeshHalo(Phi);
for (int k=1; k<Nz-1; k++){ for (int k=1; k<Nz-1; k++){
@ -544,6 +544,17 @@ void SubPhase::Full(){
gwc.K=sumReduce( Dm->Comm, wc.K); gwc.K=sumReduce( Dm->Comm, wc.K);
gwc.p=sumReduce( Dm->Comm, wc.p); gwc.p=sumReduce( Dm->Comm, wc.p);
giwn.Mn=sumReduce( Dm->Comm, iwn.Mn);
giwn.Pnx=sumReduce( Dm->Comm, iwn.Pnx);
giwn.Pny=sumReduce( Dm->Comm, iwn.Pny);
giwn.Pnz=sumReduce( Dm->Comm, iwn.Pnz);
giwn.Kn=sumReduce( Dm->Comm, iwn.Kn);
giwn.Mw=sumReduce( Dm->Comm, iwn.Mw);
giwn.Pwx=sumReduce( Dm->Comm, iwn.Pwx);
giwn.Pwy=sumReduce( Dm->Comm, iwn.Pwy);
giwn.Pwz=sumReduce( Dm->Comm, iwn.Pwz);
giwn.Kw=sumReduce( Dm->Comm, iwn.Kw);
// pressure averaging // pressure averaging
if (vol_wc_bulk > 0.0) if (vol_wc_bulk > 0.0)
wc.p = wc.p /vol_wc_bulk; wc.p = wc.p /vol_wc_bulk;

View File

@ -14,13 +14,13 @@ using namespace std;
inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
// initialize a bubble // initialize a bubble
int i,j,k,n; int i,j,k,n;
int rank = ColorModel.Mask->rank(); int rank = ColorModel.Dm->rank();
int nprocx = ColorModel.Mask->nprocx(); int nprocx = ColorModel.Dm->nprocx();
int nprocy = ColorModel.Mask->nprocy(); int nprocy = ColorModel.Dm->nprocy();
int nprocz = ColorModel.Mask->nprocz(); int nprocz = ColorModel.Dm->nprocz();
int Nx = ColorModel.Mask->Nx; int Nx = ColorModel.Dm->Nx;
int Ny = ColorModel.Mask->Ny; int Ny = ColorModel.Dm->Ny;
int Nz = ColorModel.Mask->Nz; int Nz = ColorModel.Dm->Nz;
if (rank == 0) cout << "Setting up bubble..." << endl; if (rank == 0) cout << "Setting up bubble..." << endl;
for (k=0;k<Nz;k++){ for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){ for (j=0;j<Ny;j++){
@ -34,14 +34,12 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
for (k=0;k<Nz;k++){ for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){ for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){ for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i; n = k*Nx*Ny + j*Nx + i;
int iglobal= i+(Nx-2)*ColorModel.Mask->iproc(); double iglobal= double(i+(Nx-2)*ColorModel.Dm->iproc())-double((Nx-2)*nprocx)*0.5;
int jglobal= j+(Ny-2)*ColorModel.Mask->jproc(); double jglobal= double(j+(Ny-2)*ColorModel.Dm->jproc())-double((Ny-2)*nprocy)*0.5;
int kglobal= k+(Nz-2)*ColorModel.Mask->kproc(); double kglobal= double(k+(Nz-2)*ColorModel.Dm->kproc())-double((Nz-2)*nprocz)*0.5;
// Initialize phase position field for parallel bubble test // Initialize phase position field for parallel bubble test
if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx) if ((iglobal*iglobal)+(jglobal*jglobal)+(kglobal*kglobal) < BubbleRadius*BubbleRadius){
+(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.Mask->id[n] = 2; ColorModel.Mask->id[n] = 2;
ColorModel.Mask->id[n] = 2; ColorModel.Mask->id[n] = 2;
} }
@ -50,9 +48,17 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
ColorModel.Mask->id[n]=1; ColorModel.Mask->id[n]=1;
} }
ColorModel.id[n] = ColorModel.Mask->id[n]; ColorModel.id[n] = ColorModel.Mask->id[n];
ColorModel.Dm->id[n] = ColorModel.Mask->id[n];
} }
} }
} }
FILE *OUTFILE;
char LocalRankFilename[40];
sprintf(LocalRankFilename,"Bubble.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(ColorModel.id,1,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
// initialize the phase indicator field // initialize the phase indicator field
} }
//*************************************************************************************** //***************************************************************************************