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

debugging flux

Browse Source
This commit is contained in:
JamesEMcclure 2021-06-10 10:25:08 -04:00
parent cc1a280edc
commit 9d2c46ee31
1 changed files with 54 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
models/ColorModel.cpp
View File

@ -1975,6 +1975,13 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double));
ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double));
/* DEBUG STRUCTURES */
int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz;
int N = Nx*Ny*Nz;
double * DebugMassA, *DebugMassB;
DebugMassA = new double[Np];
DebugMassB = new double[Np];
/* compute the total momentum */
vax = vay = vaz = 0.0;
vbx = vby = vbz = 0.0;
@ -2043,7 +2050,6 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
double local_momentum = sqrt(vx*vx+vy*vy+vz*vz);
if (phi > 0.0){
LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A;
if (fabs(LOCAL_MASS_CHANGE) > 0.1) LOCAL_MASS_CHANGE *= 0.1/fabs(LOCAL_MASS_CHANGE);
Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE;
Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
@ -2051,10 +2057,10 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
}
else{
LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
if (fabs(LOCAL_MASS_CHANGE) > 0.1) LOCAL_MASS_CHANGE *= 0.1/fabs(LOCAL_MASS_CHANGE);
Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE;
Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
@ -2062,6 +2068,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
DebugMassB[n] = LOCAL_MASS_CHANGE;
}
}
@ -2080,6 +2087,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE;
DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE;
}
else{
LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B;
@ -2090,11 +2098,54 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE;
DebugMassB[n] = LOCAL_MASS_CHANGE;
}
}
if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global);
/* Print out debugging info with mass update */
// initialize the array
double value;
char LocalRankFilename[40];
DoubleArray regdata(Nx,Ny,Nz);
regdata.fill(0.f);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
int idx=M.Map(i,j,k);
if (!(idx<0)){
value=DebugMassA[idx];
regdata(i,j,k)=value;
}
}
}
}
FILE *AFILE;
sprintf(LocalRankFilename,"dA.%05i.raw",M.rank);
AFILE = fopen(LocalRankFilename,"wb");
fwrite(regdata.data(),8,N,AFILE);
fclose(AFILE);
regdata.fill(0.f);
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
int idx=M.Map(i,j,k);
if (!(idx<0)){
value=DebugMassB[idx];
regdata(i,j,k)=value;
}
}
}
}
FILE *BFILE;
sprintf(LocalRankFilename,"dB.%05i.raw",M.rank);
BFILE = fopen(LocalRankFilename,"wb");
fwrite(regdata.data(),8,N,BFILE);
fclose(BFILE);
// Need to initialize Aq, Bq, Den, Phi directly
//ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));