955 lines
36 KiB
C++
955 lines
36 KiB
C++
/*
|
|
Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University
|
|
Copyright Equnior ASA
|
|
|
|
This file is part of the Open Porous Media project (OPM).
|
|
OPM is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
OPM is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
You should have received a copy of the GNU General Public License
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
#include "analysis/SubPhase.h"
|
|
|
|
// Constructor
|
|
SubPhase::SubPhase(std::shared_ptr<Domain> dm) : Dm(dm) {
|
|
Nx = dm->Nx;
|
|
Ny = dm->Ny;
|
|
Nz = dm->Nz;
|
|
Volume = (Nx - 2) * (Ny - 2) * (Nz - 2) * Dm->nprocx() * Dm->nprocy() *
|
|
Dm->nprocz() * 1.0;
|
|
|
|
morph_w = std::shared_ptr<Minkowski>(new Minkowski(Dm));
|
|
morph_n = std::shared_ptr<Minkowski>(new Minkowski(Dm));
|
|
morph_i = std::shared_ptr<Minkowski>(new Minkowski(Dm));
|
|
|
|
// Global arrays
|
|
PhaseID.resize(Nx, Ny, Nz);
|
|
PhaseID.fill(0);
|
|
Label_WP.resize(Nx, Ny, Nz);
|
|
Label_WP.fill(0);
|
|
Label_NWP.resize(Nx, Ny, Nz);
|
|
Label_NWP.fill(0);
|
|
Rho_n.resize(Nx, Ny, Nz);
|
|
Rho_n.fill(0);
|
|
Rho_w.resize(Nx, Ny, Nz);
|
|
Rho_w.fill(0);
|
|
Pressure.resize(Nx, Ny, Nz);
|
|
Pressure.fill(0);
|
|
Phi.resize(Nx, Ny, Nz);
|
|
Phi.fill(0);
|
|
DelPhi.resize(Nx, Ny, Nz);
|
|
DelPhi.fill(0);
|
|
Vel_x.resize(Nx, Ny, Nz);
|
|
Vel_x.fill(0); // Gradient of the phase indicator field
|
|
Vel_y.resize(Nx, Ny, Nz);
|
|
Vel_y.fill(0);
|
|
Vel_z.resize(Nx, Ny, Nz);
|
|
Vel_z.fill(0);
|
|
Dissipation.resize(Nx, Ny, Nz);
|
|
Dissipation.fill(0);
|
|
SDs.resize(Nx, Ny, Nz);
|
|
SDs.fill(0);
|
|
//.........................................
|
|
|
|
//.........................................
|
|
if (Dm->rank() == 0) {
|
|
bool WriteHeader = false;
|
|
SUBPHASE = fopen("subphase.csv", "r");
|
|
if (SUBPHASE != NULL)
|
|
fclose(SUBPHASE);
|
|
else
|
|
WriteHeader = true;
|
|
|
|
SUBPHASE = fopen("subphase.csv", "a+");
|
|
if (WriteHeader) {
|
|
// If timelog is empty, write a short header to list the averages
|
|
fprintf(SUBPHASE, "time rn rw nun nuw Fx Fy Fz iftwn wet ");
|
|
fprintf(SUBPHASE, "pwc pwd pnc pnd "); // pressures
|
|
fprintf(SUBPHASE, "Mwc Mwd Mwi Mnc Mnd Mni Msw Msn "); // mass
|
|
fprintf(
|
|
SUBPHASE,
|
|
"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x "); // momentum
|
|
fprintf(SUBPHASE,
|
|
"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y ");
|
|
fprintf(SUBPHASE,
|
|
"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z ");
|
|
fprintf(SUBPHASE, "Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
|
|
fprintf(SUBPHASE, "Dwc Dwd Dnc Dnd "); // viscous dissipation
|
|
fprintf(SUBPHASE, "Vwc Awc Hwc Xwc "); // wc region
|
|
fprintf(SUBPHASE, "Vwd Awd Hwd Xwd Nwd "); // wd region
|
|
fprintf(SUBPHASE, "Vnc Anc Hnc Xnc "); // nc region
|
|
fprintf(SUBPHASE, "Vnd And Hnd Xnd Nnd "); // nd regionin
|
|
fprintf(SUBPHASE, "Vi Ai Hi Xi "); // interface region
|
|
fprintf(SUBPHASE, "Vic Aic Hic Xic Nic\n"); // interface region
|
|
|
|
// stress tensor?
|
|
}
|
|
|
|
} else {
|
|
char LocalRankString[8];
|
|
sprintf(LocalRankString, "%05d", Dm->rank());
|
|
char LocalRankFilename[40];
|
|
sprintf(LocalRankFilename, "%s%s", "subphase.csv.", LocalRankString);
|
|
SUBPHASE = fopen(LocalRankFilename, "a+");
|
|
|
|
fprintf(SUBPHASE, "time rn rw nun nuw Fx Fy Fz iftwn wet ");
|
|
fprintf(SUBPHASE, "pwc pwd pnc pnd "); // pressures
|
|
fprintf(SUBPHASE, "Mwc Mwd Mwi Mnc Mnd Mni Msw Msn "); // mass
|
|
fprintf(SUBPHASE,
|
|
"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x "); // momentum
|
|
fprintf(SUBPHASE, "Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y ");
|
|
fprintf(SUBPHASE, "Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z ");
|
|
fprintf(SUBPHASE, "Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
|
|
fprintf(SUBPHASE, "Dwc Dwd Dnc Dnd "); // viscous dissipation
|
|
fprintf(SUBPHASE, "Vwc Awc Hwc Xwc "); // wc region
|
|
fprintf(SUBPHASE, "Vwd Awd Hwd Xwd Nwd "); // wd region
|
|
fprintf(SUBPHASE, "Vnc Anc Hnc Xnc "); // nc region
|
|
fprintf(SUBPHASE, "Vnd And Hnd Xnd Nnd "); // nd region
|
|
fprintf(SUBPHASE, "Vi Ai Hi Xi "); // interface region
|
|
fprintf(SUBPHASE, "Vic Aic Hic Xic Nic\n"); // interface region
|
|
}
|
|
|
|
if (Dm->rank() == 0) {
|
|
bool WriteHeader = false;
|
|
TIMELOG = fopen("timelog.csv", "r");
|
|
if (TIMELOG != NULL)
|
|
fclose(TIMELOG);
|
|
else
|
|
WriteHeader = true;
|
|
|
|
TIMELOG = fopen("timelog.csv", "a+");
|
|
if (WriteHeader) {
|
|
// If timelog is empty, write a short header to list the averages
|
|
fprintf(TIMELOG,
|
|
"sw krw krn krwf krnf vw vn force pw pn wet peff\n");
|
|
}
|
|
}
|
|
}
|
|
|
|
// Destructor
|
|
SubPhase::~SubPhase() {
|
|
if (SUBPHASE != NULL) {
|
|
fclose(SUBPHASE);
|
|
}
|
|
}
|
|
|
|
void SubPhase::Write(int timestep) {
|
|
if (Dm->rank() == 0) {
|
|
fprintf(SUBPHASE, "%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",
|
|
timestep, rho_n, rho_w, nu_n, nu_w, Fx, Fy, Fz, gamma_wn,
|
|
total_wetting_interaction_global);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", gwc.p, gwd.p, gnc.p, gnd.p);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", gwc.M,
|
|
gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn, gifs.Mw, gifs.Mn);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", gwc.Px,
|
|
gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx, gifs.Pwx, gifs.Pnx);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", gwc.Py,
|
|
gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny, gifs.Pwy, gifs.Pny);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", gwc.Pz,
|
|
gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz, gifs.Pwz, gifs.Pnz);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g ", gwc.K, gwd.K,
|
|
giwn.Kw, gnc.K, gnd.K, giwn.Kn);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", gwc.visc, gwd.visc, gnc.visc,
|
|
gnd.visc);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", gwc.V, gwc.A, gwc.H, gwc.X);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %i ", gwd.V, gwd.A, gwd.H, gwd.X,
|
|
gwd.Nc);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", gnc.V, gnc.A, gnc.H, gnc.X);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %i ", gnd.V, gnd.A, gnd.H, gnd.X,
|
|
gnd.Nc);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", giwn.V, giwn.A, giwn.H,
|
|
giwn.X);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %i\n", giwnc.V, giwnc.A, giwnc.H,
|
|
giwnc.X, giwnc.Nc);
|
|
fflush(SUBPHASE);
|
|
} else {
|
|
fprintf(SUBPHASE, "%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",
|
|
timestep, rho_n, rho_w, nu_n, nu_w, Fx, Fy, Fz, gamma_wn,
|
|
total_wetting_interaction);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", wc.p, wd.p, nc.p, nd.p);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", wc.M,
|
|
wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn, ifs.Mw, ifs.Mn);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", wc.Px,
|
|
wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx, ifs.Pwx, ifs.Pnx);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", wc.Py,
|
|
wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny, ifs.Pwy, ifs.Pny);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ", wc.Pz,
|
|
wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz, ifs.Pwz, ifs.Pnz);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %.8g %.8g ", wc.K, wd.K, iwn.Kw,
|
|
nc.K, nd.K, iwn.Kn);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", wc.visc, wd.visc, nc.visc,
|
|
nd.visc);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", wc.V, wc.A, wc.H, wc.X);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %i ", wd.V, wd.A, wd.H, wd.X,
|
|
wd.Nc);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", nc.V, nc.A, nc.H, nc.X);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g %i ", nd.V, nd.A, nd.H, nd.X,
|
|
nd.Nc);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g ", iwn.V, iwn.A, iwn.H, iwn.X);
|
|
fprintf(SUBPHASE, "%.8g %.8g %.8g %.8g\n", iwnc.V, iwnc.A, iwnc.H,
|
|
iwnc.X);
|
|
}
|
|
}
|
|
|
|
void SubPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB,
|
|
double force_x, double force_y, double force_z,
|
|
double alpha, double B) {
|
|
Fx = force_x;
|
|
Fy = force_y;
|
|
Fz = force_z;
|
|
rho_n = rhoA;
|
|
rho_w = rhoB;
|
|
nu_n = (tauA - 0.5) / 3.f;
|
|
nu_w = (tauB - 0.5) / 3.f;
|
|
gamma_wn = 5.796 * alpha;
|
|
beta = B;
|
|
}
|
|
|
|
void SubPhase::Basic() {
|
|
int i, j, k, n, imin, jmin, kmin, kmax;
|
|
|
|
// If external boundary conditions are set, do not average over the inlet
|
|
kmin = 1;
|
|
kmax = Nz - 1;
|
|
imin = jmin = 1;
|
|
|
|
nb.reset();
|
|
wb.reset();
|
|
iwn.reset();
|
|
|
|
double count_w = 0.0;
|
|
double count_n = 0.0;
|
|
|
|
/* compute the laplacian */
|
|
Dm->CommunicateMeshHalo(Phi);
|
|
for (int k = 1; k < Nz - 1; k++) {
|
|
for (int j = 1; j < Ny - 1; j++) {
|
|
for (int i = 1; i < Nx - 1; i++) {
|
|
// Compute all of the derivatives using finite differences
|
|
double fx = 0.5 * (Phi(i + 1, j, k) - Phi(i - 1, j, k));
|
|
double fy = 0.5 * (Phi(i, j + 1, k) - Phi(i, j - 1, k));
|
|
double fz = 0.5 * (Phi(i, j, k + 1) - Phi(i, j, k - 1));
|
|
DelPhi(i, j, k) = sqrt(fx * fx + fy * fy + fz * fz);
|
|
}
|
|
}
|
|
}
|
|
Dm->CommunicateMeshHalo(DelPhi);
|
|
|
|
for (k = 0; k < Nz; k++) {
|
|
for (j = 0; j < Ny; j++) {
|
|
for (i = 0; i < Nx; i++) {
|
|
n = k * Nx * Ny + j * Nx + i;
|
|
// Compute volume averages
|
|
if (Dm->id[n] > 0) {
|
|
// compute density
|
|
double nA = Rho_n(n);
|
|
double nB = Rho_w(n);
|
|
double phi = (nA - nB) / (nA + nB);
|
|
Phi(n) = phi;
|
|
}
|
|
if (Phi(n) != Phi(n)) {
|
|
// check for NaN
|
|
Phi(n) = 0.0;
|
|
//printf("Nan at %i %i %i \n",i,j,k);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for (k = kmin; k < kmax; k++) {
|
|
for (j = jmin; j < Ny - 1; j++) {
|
|
for (i = imin; i < Nx - 1; i++) {
|
|
n = k * Nx * Ny + j * Nx + i;
|
|
// Compute volume averages
|
|
if (Dm->id[n] > 0) {
|
|
// compute density
|
|
double nA = Rho_n(n);
|
|
double nB = Rho_w(n);
|
|
double phi = (nA - nB) / (nA + nB);
|
|
if (phi > 0.0) {
|
|
nA = 1.0;
|
|
nb.V += 1.0;
|
|
nb.M += nA * rho_n;
|
|
// velocity
|
|
nb.Px += rho_n * nA * Vel_x(n);
|
|
nb.Py += rho_n * nA * Vel_y(n);
|
|
nb.Pz += rho_n * nA * Vel_z(n);
|
|
} else {
|
|
nB = 1.0;
|
|
wb.M += nB * rho_w;
|
|
wb.V += 1.0;
|
|
|
|
// velocity
|
|
wb.Px += rho_w * nB * Vel_x(n);
|
|
wb.Py += rho_w * nB * Vel_y(n);
|
|
wb.Pz += rho_w * nB * Vel_z(n);
|
|
}
|
|
if (phi > 0.99) {
|
|
nb.p += Pressure(n);
|
|
count_n += 1.0;
|
|
} else if (phi < -0.99) {
|
|
wb.p += Pressure(n);
|
|
count_w += 1.0;
|
|
}
|
|
/* compute the film contribution */
|
|
else if (SDs(i, j, k) < 2.0) {
|
|
if (phi > 0.0) {
|
|
nA = 1.0;
|
|
iwn.V += 1.0;
|
|
iwn.Mn += nA * rho_n;
|
|
// velocity
|
|
iwn.Pnx += rho_n * nA * Vel_x(n);
|
|
iwn.Pny += rho_n * nA * Vel_y(n);
|
|
iwn.Pnz += rho_n * nA * Vel_z(n);
|
|
} else {
|
|
nB = 1.0;
|
|
iwn.Mw += nB * rho_w;
|
|
iwn.V += 1.0;
|
|
|
|
iwn.Pwx += rho_w * nB * Vel_x(n);
|
|
iwn.Pwy += rho_w * nB * Vel_y(n);
|
|
iwn.Pwz += rho_w * nB * Vel_z(n);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
total_wetting_interaction = count_wetting_interaction = 0.0;
|
|
total_wetting_interaction_global = count_wetting_interaction_global = 0.0;
|
|
for (k = kmin; k < kmax; k++) {
|
|
for (j = jmin; j < Ny - 1; j++) {
|
|
for (i = imin; i < Nx - 1; i++) {
|
|
n = k * Nx * Ny + j * Nx + i;
|
|
// compute contribution of wetting terms (within two voxels of solid)
|
|
if (Dm->id[n] > 0 && SDs(i, j, k) < 2.0) {
|
|
count_wetting_interaction += 1.0;
|
|
total_wetting_interaction += DelPhi(i, j, k);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
total_wetting_interaction_global =
|
|
Dm->Comm.sumReduce(total_wetting_interaction);
|
|
count_wetting_interaction_global =
|
|
Dm->Comm.sumReduce(count_wetting_interaction);
|
|
|
|
gwb.V = Dm->Comm.sumReduce(wb.V);
|
|
gnb.V = Dm->Comm.sumReduce(nb.V);
|
|
gwb.M = Dm->Comm.sumReduce(wb.M);
|
|
gnb.M = Dm->Comm.sumReduce(nb.M);
|
|
gwb.Px = Dm->Comm.sumReduce(wb.Px);
|
|
gwb.Py = Dm->Comm.sumReduce(wb.Py);
|
|
gwb.Pz = Dm->Comm.sumReduce(wb.Pz);
|
|
gnb.Px = Dm->Comm.sumReduce(nb.Px);
|
|
gnb.Py = Dm->Comm.sumReduce(nb.Py);
|
|
gnb.Pz = Dm->Comm.sumReduce(nb.Pz);
|
|
|
|
giwn.Mw = Dm->Comm.sumReduce(iwn.Mw);
|
|
giwn.Pwx = Dm->Comm.sumReduce(iwn.Pwx);
|
|
giwn.Pwy = Dm->Comm.sumReduce(iwn.Pwy);
|
|
giwn.Pwz = Dm->Comm.sumReduce(iwn.Pwz);
|
|
|
|
giwn.Mn = Dm->Comm.sumReduce(iwn.Mn);
|
|
giwn.Pnx = Dm->Comm.sumReduce(iwn.Pnx);
|
|
giwn.Pny = Dm->Comm.sumReduce(iwn.Pny);
|
|
giwn.Pnz = Dm->Comm.sumReduce(iwn.Pnz);
|
|
|
|
count_w = Dm->Comm.sumReduce(count_w);
|
|
count_n = Dm->Comm.sumReduce(count_n);
|
|
if (count_w > 0.0)
|
|
gwb.p = Dm->Comm.sumReduce(wb.p) / count_w;
|
|
else
|
|
gwb.p = 0.0;
|
|
if (count_n > 0.0)
|
|
gnb.p = Dm->Comm.sumReduce(nb.p) / count_n;
|
|
else
|
|
gnb.p = 0.0;
|
|
|
|
// check for NaN
|
|
bool err = false;
|
|
if (gwb.V != gwb.V)
|
|
err = true;
|
|
if (gnb.V != gnb.V)
|
|
err = true;
|
|
if (gwb.p != gwb.p)
|
|
err = true;
|
|
if (gnb.p != gnb.p)
|
|
err = true;
|
|
if (gwb.Px != gwb.Px)
|
|
err = true;
|
|
if (gwb.Py != gwb.Py)
|
|
err = true;
|
|
if (gwb.Pz != gwb.Pz)
|
|
err = true;
|
|
if (gnb.Px != gnb.Px)
|
|
err = true;
|
|
if (gnb.Py != gnb.Py)
|
|
err = true;
|
|
if (gnb.Pz != gnb.Pz)
|
|
err = true;
|
|
|
|
if (Dm->rank() == 0) {
|
|
/* align flow direction based on total mass flux */
|
|
double dir_x = gwb.Px + gnb.Px;
|
|
double dir_y = gwb.Py + gnb.Py;
|
|
double dir_z = gwb.Pz + gnb.Pz;
|
|
double flow_magnitude = dir_x * dir_x + dir_y * dir_y + dir_z * dir_z;
|
|
double force_mag = sqrt(Fx * Fx + Fy * Fy + Fz * Fz);
|
|
if (force_mag > 0.0) {
|
|
dir_x = Fx / force_mag;
|
|
dir_y = Fy / force_mag;
|
|
dir_z = Fz / force_mag;
|
|
} else {
|
|
dir_x /= flow_magnitude;
|
|
dir_y /= flow_magnitude;
|
|
dir_z /= flow_magnitude;
|
|
}
|
|
if (Dm->BoundaryCondition == 1 || Dm->BoundaryCondition == 2 ||
|
|
Dm->BoundaryCondition == 3 || Dm->BoundaryCondition == 4) {
|
|
// compute the pressure drop
|
|
double pressure_drop = (Pressure(Nx * Ny + Nx + 1) - 1.0) / 3.0;
|
|
double length = ((Nz - 2) * Dm->nprocz());
|
|
force_mag -= pressure_drop / length;
|
|
}
|
|
if (force_mag == 0.0 && flow_magnitude == 0.0) {
|
|
// default to z direction
|
|
dir_x = 0.0;
|
|
dir_y = 0.0;
|
|
dir_z = 1.0;
|
|
force_mag = 1.0;
|
|
}
|
|
double Porosity = (gwb.V + gnb.V) / Dm->Volume;
|
|
double saturation = gwb.V / (gwb.V + gnb.V);
|
|
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;
|
|
|
|
/* contribution from water films */
|
|
double water_film_flow_rate =
|
|
gwb.V * (giwn.Pwx * dir_x + giwn.Pwy * dir_y + giwn.Pwz * dir_z) /
|
|
gwb.M / Dm->Volume;
|
|
double not_water_film_flow_rate =
|
|
gnb.V * (giwn.Pnx * dir_x + giwn.Pny * dir_y + giwn.Pnz * 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;
|
|
double h = Dm->voxel_length;
|
|
double krn = h * h * nu_n * Porosity * not_water_flow_rate / force_mag;
|
|
double krw = h * h * nu_w * Porosity * water_flow_rate / force_mag;
|
|
/* not counting films */
|
|
double krnf = krn - h * h * nu_n * Porosity * not_water_film_flow_rate /
|
|
force_mag;
|
|
double krwf =
|
|
krw - h * h * nu_w * Porosity * water_film_flow_rate / force_mag;
|
|
double eff_pressure = 1.0 / (krn + krw); // effective pressure drop
|
|
|
|
fprintf(TIMELOG,
|
|
"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",
|
|
saturation, krw, krn, krwf, krnf, h * water_flow_rate,
|
|
h * not_water_flow_rate, force_mag, gwb.p, gnb.p,
|
|
total_wetting_interaction_global, eff_pressure);
|
|
fflush(TIMELOG);
|
|
}
|
|
if (err == true) {
|
|
// exception if simulation produceds NaN
|
|
printf("SubPhase.cpp: NaN encountered, may need to check simulation "
|
|
"parameters \n");
|
|
}
|
|
ASSERT(err == false);
|
|
}
|
|
|
|
inline void InterfaceTransportMeasures(double beta, double rA, double rB,
|
|
double nA, double nB, double nx,
|
|
double ny, double nz, double ux,
|
|
double uy, double uz, interface &I) {
|
|
|
|
double A1, A2, A3, A4, A5, A6;
|
|
double B1, B2, B3, B4, B5, B6;
|
|
double nAB, delta;
|
|
double phi = (nA - nB) / (nA + nB);
|
|
// Instantiate mass transport distributions
|
|
// Stationary value - distribution 0
|
|
nAB = 1.0 / (nA + nB);
|
|
//...............................................
|
|
// q = 0,2,4
|
|
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
|
|
delta = beta * nA * nB * nAB * 0.1111111111111111 * nx;
|
|
if (!(nA * nB * nAB > 0))
|
|
delta = 0;
|
|
A1 = nA * (0.1111111111111111 * (1 + 4.5 * ux)) + delta;
|
|
B1 = nB * (0.1111111111111111 * (1 + 4.5 * ux)) - delta;
|
|
A2 = nA * (0.1111111111111111 * (1 - 4.5 * ux)) - delta;
|
|
B2 = nB * (0.1111111111111111 * (1 - 4.5 * ux)) + delta;
|
|
|
|
//...............................................
|
|
// Cq = {0,1,0}
|
|
delta = beta * nA * nB * nAB * 0.1111111111111111 * ny;
|
|
if (!(nA * nB * nAB > 0))
|
|
delta = 0;
|
|
A3 = nA * (0.1111111111111111 * (1 + 4.5 * uy)) + delta;
|
|
B3 = nB * (0.1111111111111111 * (1 + 4.5 * uy)) - delta;
|
|
A4 = nA * (0.1111111111111111 * (1 - 4.5 * uy)) - delta;
|
|
B4 = nB * (0.1111111111111111 * (1 - 4.5 * uy)) + delta;
|
|
|
|
//...............................................
|
|
// q = 4
|
|
// Cq = {0,0,1}
|
|
delta = beta * nA * nB * nAB * 0.1111111111111111 * nz;
|
|
if (!(nA * nB * nAB > 0))
|
|
delta = 0;
|
|
A5 = nA * (0.1111111111111111 * (1 + 4.5 * uz)) + delta;
|
|
B5 = nB * (0.1111111111111111 * (1 + 4.5 * uz)) - delta;
|
|
A6 = nA * (0.1111111111111111 * (1 - 4.5 * uz)) - delta;
|
|
B6 = nB * (0.1111111111111111 * (1 - 4.5 * uz)) + delta;
|
|
|
|
double unx = (A1 - A2);
|
|
double uny = (A3 - A4);
|
|
double unz = (A5 - A6);
|
|
double uwx = (B1 - B2);
|
|
double uwy = (B3 - B4);
|
|
double uwz = (B5 - B6);
|
|
/*
|
|
I.Mn += rA*nA;
|
|
I.Mw += rB*nB;
|
|
I.Pnx += rA*nA*unx;
|
|
I.Pny += rA*nA*uny;
|
|
I.Pnz += rA*nA*unz;
|
|
I.Pwx += rB*nB*uwx;
|
|
I.Pwy += rB*nB*uwy;
|
|
I.Pwz += rB*nB*uwz;
|
|
I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz);
|
|
I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz);
|
|
*/
|
|
if (phi > 0.0) {
|
|
I.Mn += rA;
|
|
I.Pnx += rA * ux;
|
|
I.Pny += rA * uy;
|
|
I.Pnz += rA * uz;
|
|
} else {
|
|
I.Mw += rB;
|
|
I.Pwx += rB * ux;
|
|
I.Pwy += rB * uy;
|
|
I.Pwz += rB * uz;
|
|
}
|
|
I.Kn += rA * nA * (unx * unx + uny * uny + unz * unz);
|
|
I.Kw += rB * nB * (uwx * uwx + uwy * uwy + uwz * uwz);
|
|
}
|
|
|
|
void SubPhase::Full() {
|
|
int i, j, k, n, imin, jmin, kmin, kmax;
|
|
|
|
// If external boundary conditions are set, do not average over the inlet
|
|
kmin = 1;
|
|
kmax = Nz - 1;
|
|
imin = jmin = 1;
|
|
|
|
nd.reset();
|
|
nc.reset();
|
|
wd.reset();
|
|
wc.reset();
|
|
iwn.reset();
|
|
iwnc.reset();
|
|
ifs.reset();
|
|
|
|
Dm->CommunicateMeshHalo(Phi);
|
|
for (int k = 1; k < Nz - 1; k++) {
|
|
for (int j = 1; j < Ny - 1; j++) {
|
|
for (int i = 1; i < Nx - 1; i++) {
|
|
// Compute all of the derivatives using finite differences
|
|
double fx = 0.5 * (Phi(i + 1, j, k) - Phi(i - 1, j, k));
|
|
double fy = 0.5 * (Phi(i, j + 1, k) - Phi(i, j - 1, k));
|
|
double fz = 0.5 * (Phi(i, j, k + 1) - Phi(i, j, k - 1));
|
|
DelPhi(i, j, k) = sqrt(fx * fx + fy * fy + fz * fz);
|
|
}
|
|
}
|
|
}
|
|
Dm->CommunicateMeshHalo(DelPhi);
|
|
|
|
Dm->CommunicateMeshHalo(Vel_x);
|
|
Dm->CommunicateMeshHalo(Vel_y);
|
|
Dm->CommunicateMeshHalo(Vel_z);
|
|
for (int k = 1; k < Nz - 1; k++) {
|
|
for (int j = 1; j < Ny - 1; j++) {
|
|
for (int i = 1; i < Nx - 1; i++) {
|
|
// Compute velocity gradients using finite differences
|
|
double phi = Phi(i, j, k);
|
|
double nu = nu_n + 0.5 * (1.0 - phi) * (nu_w - nu_n);
|
|
double rho = rho_n + 0.5 * (1.0 - phi) * (rho_w - rho_n);
|
|
double ux = 0.5 * (Vel_x(i + 1, j, k) - Vel_x(i - 1, j, k));
|
|
double uy = 0.5 * (Vel_x(i, j + 1, k) - Vel_x(i, j - 1, k));
|
|
double uz = 0.5 * (Vel_x(i, j, k + 1) - Vel_x(i, j, k - 1));
|
|
double vx = 0.5 * (Vel_y(i + 1, j, k) - Vel_y(i - 1, j, k));
|
|
double vy = 0.5 * (Vel_y(i, j + 1, k) - Vel_y(i, j - 1, k));
|
|
double vz = 0.5 * (Vel_y(i, j, k + 1) - Vel_y(i, j, k - 1));
|
|
double wx = 0.5 * (Vel_z(i + 1, j, k) - Vel_z(i - 1, j, k));
|
|
double wy = 0.5 * (Vel_z(i, j + 1, k) - Vel_z(i, j - 1, k));
|
|
double wz = 0.5 * (Vel_z(i, j, k + 1) - Vel_z(i, j, k - 1));
|
|
if (SDs(i, j, k) > 2.0) {
|
|
Dissipation(i, j, k) = 2 * rho * nu *
|
|
(ux * ux + vy * vy + wz * wz +
|
|
0.5 * (vx + uy) * (vx + uy) +
|
|
0.5 * (vz + wy) * (vz + wy) +
|
|
0.5 * (uz + wx) * (uz + wx));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Set up geometric analysis of each region */
|
|
|
|
// non-wetting
|
|
for (k = 0; k < Nz; k++) {
|
|
for (j = 0; j < Ny; j++) {
|
|
for (i = 0; i < Nx; i++) {
|
|
n = k * Nx * Ny + j * Nx + i;
|
|
if (SDs(n) <= 0.0) {
|
|
// Solid phase
|
|
morph_n->id(i, j, k) = 1;
|
|
|
|
} else if (Phi(n) > 0.0) {
|
|
// non-wetting phase
|
|
morph_n->id(i, j, k) = 0;
|
|
} else {
|
|
// wetting phase
|
|
morph_n->id(i, j, k) = 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// measure the whole object
|
|
morph_n->MeasureObject(); //0.5/beta,Phi);
|
|
nd.V = morph_n->V();
|
|
nd.A = morph_n->A();
|
|
nd.H = morph_n->H();
|
|
nd.X = morph_n->X();
|
|
// measure only the connected part
|
|
nd.Nc = morph_n->MeasureConnectedPathway(); //0.5/beta,Phi);
|
|
nc.V = morph_n->V();
|
|
nc.A = morph_n->A();
|
|
nc.H = morph_n->H();
|
|
nc.X = morph_n->X();
|
|
// update disconnected part
|
|
nd.V -= nc.V;
|
|
nd.A -= nc.A;
|
|
nd.H -= nc.H;
|
|
nd.X -= nc.X;
|
|
|
|
// compute global entities
|
|
gnc.V = Dm->Comm.sumReduce(nc.V);
|
|
gnc.A = Dm->Comm.sumReduce(nc.A);
|
|
gnc.H = Dm->Comm.sumReduce(nc.H);
|
|
gnc.X = Dm->Comm.sumReduce(nc.X);
|
|
gnd.V = Dm->Comm.sumReduce(nd.V);
|
|
gnd.A = Dm->Comm.sumReduce(nd.A);
|
|
gnd.H = Dm->Comm.sumReduce(nd.H);
|
|
gnd.X = Dm->Comm.sumReduce(nd.X);
|
|
gnd.Nc = nd.Nc;
|
|
// wetting
|
|
for (k = 0; k < Nz; k++) {
|
|
for (j = 0; j < Ny; j++) {
|
|
for (i = 0; i < Nx; i++) {
|
|
n = k * Nx * Ny + j * Nx + i;
|
|
if (SDs(n) <= 0.0) {
|
|
// Solid phase
|
|
morph_w->id(i, j, k) = 1;
|
|
|
|
} else if (Phi(n) < 0.0) {
|
|
// wetting phase
|
|
morph_w->id(i, j, k) = 0;
|
|
} else {
|
|
// non-wetting phase
|
|
morph_w->id(i, j, k) = 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
morph_w->MeasureObject(); //-0.5/beta,Phi);
|
|
wd.V = morph_w->V();
|
|
wd.A = morph_w->A();
|
|
wd.H = morph_w->H();
|
|
wd.X = morph_w->X();
|
|
// measure only the connected part
|
|
wd.Nc = morph_w->MeasureConnectedPathway(); //-0.5/beta,Phi);
|
|
wc.V = morph_w->V();
|
|
wc.A = morph_w->A();
|
|
wc.H = morph_w->H();
|
|
wc.X = morph_w->X();
|
|
// update disconnected part
|
|
wd.V -= wc.V;
|
|
wd.A -= wc.A;
|
|
wd.H -= wc.H;
|
|
wd.X -= wc.X;
|
|
// compute global entities
|
|
gwc.V = Dm->Comm.sumReduce(wc.V);
|
|
gwc.A = Dm->Comm.sumReduce(wc.A);
|
|
gwc.H = Dm->Comm.sumReduce(wc.H);
|
|
gwc.X = Dm->Comm.sumReduce(wc.X);
|
|
gwd.V = Dm->Comm.sumReduce(wd.V);
|
|
gwd.A = Dm->Comm.sumReduce(wd.A);
|
|
gwd.H = Dm->Comm.sumReduce(wd.H);
|
|
gwd.X = Dm->Comm.sumReduce(wd.X);
|
|
gwd.Nc = wd.Nc;
|
|
|
|
/* Set up geometric analysis of interface region */
|
|
for (k = 0; k < Nz; k++) {
|
|
for (j = 0; j < Ny; j++) {
|
|
for (i = 0; i < Nx; i++) {
|
|
n = k * Nx * Ny + j * Nx + i;
|
|
if (SDs(n) <= 0.0) {
|
|
// Solid phase
|
|
morph_i->id(i, j, k) = 1;
|
|
} else if (DelPhi(n) > 1e-4) {
|
|
// interface
|
|
morph_i->id(i, j, k) = 0;
|
|
} else {
|
|
// not interface
|
|
morph_i->id(i, j, k) = 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
morph_i->MeasureObject();
|
|
iwn.V = morph_i->V();
|
|
iwn.A = morph_i->A();
|
|
iwn.H = morph_i->H();
|
|
iwn.X = morph_i->X();
|
|
giwn.V = Dm->Comm.sumReduce(iwn.V);
|
|
giwn.A = Dm->Comm.sumReduce(iwn.A);
|
|
giwn.H = Dm->Comm.sumReduce(iwn.H);
|
|
giwn.X = Dm->Comm.sumReduce(iwn.X);
|
|
// measure only the connected part
|
|
iwnc.Nc = morph_i->MeasureConnectedPathway();
|
|
iwnc.V = morph_i->V();
|
|
iwnc.A = morph_i->A();
|
|
iwnc.H = morph_i->H();
|
|
iwnc.X = morph_i->X();
|
|
giwnc.V = Dm->Comm.sumReduce(iwnc.V);
|
|
giwnc.A = Dm->Comm.sumReduce(iwnc.A);
|
|
giwnc.H = Dm->Comm.sumReduce(iwnc.H);
|
|
giwnc.X = Dm->Comm.sumReduce(iwnc.X);
|
|
giwnc.Nc = iwnc.Nc;
|
|
|
|
double vol_nc_bulk = 0.0;
|
|
double vol_wc_bulk = 0.0;
|
|
double vol_nd_bulk = 0.0;
|
|
double vol_wd_bulk = 0.0;
|
|
for (k = kmin; k < kmax; k++) {
|
|
for (j = jmin; j < Ny - 1; j++) {
|
|
for (i = imin; i < Nx - 1; i++) {
|
|
n = k * Nx * Ny + j * Nx + i;
|
|
// Compute volume averages
|
|
if (SDs(n) > 0.0) {
|
|
// compute density
|
|
double nA = Rho_n(n);
|
|
double nB = Rho_w(n);
|
|
double phi = (nA - nB) / (nA + nB);
|
|
double ux = Vel_x(n);
|
|
double uy = Vel_y(n);
|
|
double uz = Vel_z(n);
|
|
double visc = Dissipation(n);
|
|
|
|
if (DelPhi(n) > 1e-3) {
|
|
// get the normal vector
|
|
double nx = 0.5 * (Phi(i + 1, j, k) - Phi(i - 1, j, k));
|
|
double ny = 0.5 * (Phi(i, j + 1, k) - Phi(i, j - 1, k));
|
|
double nz = 0.5 * (Phi(i, j, k + 1) - Phi(i, j, k - 1));
|
|
if (SDs(n) > 2.5) {
|
|
// not a film region
|
|
InterfaceTransportMeasures(beta, rho_w, rho_n, nA,
|
|
nB, nx, ny, nz, ux, uy,
|
|
uz, iwn);
|
|
} else {
|
|
// films that are close to the wetting fluid
|
|
if (morph_w->distance(i, j, k) < 2.5 && phi > 0.0) {
|
|
ifs.Mw += rho_w;
|
|
ifs.Pwx += rho_w * ux;
|
|
ifs.Pwy += rho_w * uy;
|
|
ifs.Pwz += rho_w * uz;
|
|
}
|
|
// films that are close to the NWP
|
|
if (morph_n->distance(i, j, k) < 2.5 && phi < 0.0) {
|
|
ifs.Mn += rho_n;
|
|
ifs.Pnx += rho_n * ux;
|
|
ifs.Pny += rho_n * uy;
|
|
ifs.Pnz += rho_n * uz;
|
|
}
|
|
}
|
|
} else if (phi > 0.0) {
|
|
if (morph_n->label(i, j, k) > 0) {
|
|
vol_nd_bulk += 1.0;
|
|
nd.p += Pressure(n);
|
|
} else {
|
|
vol_nc_bulk += 1.0;
|
|
nc.p += Pressure(n);
|
|
}
|
|
} else {
|
|
// water region
|
|
if (morph_w->label(i, j, k) > 0) {
|
|
vol_wd_bulk += 1.0;
|
|
wd.p += Pressure(n);
|
|
} else {
|
|
vol_wc_bulk += 1.0;
|
|
wc.p += Pressure(n);
|
|
}
|
|
}
|
|
if (phi > 0.0) {
|
|
if (morph_n->label(i, j, k) > 0) {
|
|
nA = 1.0;
|
|
nd.M += nA * rho_n;
|
|
nd.Px += nA * rho_n * ux;
|
|
nd.Py += nA * rho_n * uy;
|
|
nd.Pz += nA * rho_n * uz;
|
|
nd.K += nA * rho_n * (ux * ux + uy * uy + uz * uz);
|
|
nd.visc += visc;
|
|
} else {
|
|
nA = 1.0;
|
|
nc.M += nA * rho_n;
|
|
nc.Px += nA * rho_n * ux;
|
|
nc.Py += nA * rho_n * uy;
|
|
nc.Pz += nA * rho_n * uz;
|
|
nc.K += nA * rho_n * (ux * ux + uy * uy + uz * uz);
|
|
nc.visc += visc;
|
|
}
|
|
} else {
|
|
// water region
|
|
if (morph_w->label(i, j, k) > 0) {
|
|
nB = 1.0;
|
|
wd.M += nB * rho_w;
|
|
wd.Px += nB * rho_w * ux;
|
|
wd.Py += nB * rho_w * uy;
|
|
wd.Pz += nB * rho_w * uz;
|
|
wd.K += nB * rho_w * (ux * ux + uy * uy + uz * uz);
|
|
wd.visc += visc;
|
|
} else {
|
|
nB = 1.0;
|
|
wc.M += nB * rho_w;
|
|
wc.Px += nB * rho_w * ux;
|
|
wc.Py += nB * rho_w * uy;
|
|
wc.Pz += nB * rho_w * uz;
|
|
wc.K += nB * rho_w * (ux * ux + uy * uy + uz * uz);
|
|
wc.visc += visc;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
gnd.M = Dm->Comm.sumReduce(nd.M);
|
|
gnd.Px = Dm->Comm.sumReduce(nd.Px);
|
|
gnd.Py = Dm->Comm.sumReduce(nd.Py);
|
|
gnd.Pz = Dm->Comm.sumReduce(nd.Pz);
|
|
gnd.K = Dm->Comm.sumReduce(nd.K);
|
|
gnd.visc = Dm->Comm.sumReduce(nd.visc);
|
|
|
|
gwd.M = Dm->Comm.sumReduce(wd.M);
|
|
gwd.Px = Dm->Comm.sumReduce(wd.Px);
|
|
gwd.Py = Dm->Comm.sumReduce(wd.Py);
|
|
gwd.Pz = Dm->Comm.sumReduce(wd.Pz);
|
|
gwd.K = Dm->Comm.sumReduce(wd.K);
|
|
gwd.visc = Dm->Comm.sumReduce(wd.visc);
|
|
|
|
gnc.M = Dm->Comm.sumReduce(nc.M);
|
|
gnc.Px = Dm->Comm.sumReduce(nc.Px);
|
|
gnc.Py = Dm->Comm.sumReduce(nc.Py);
|
|
gnc.Pz = Dm->Comm.sumReduce(nc.Pz);
|
|
gnc.K = Dm->Comm.sumReduce(nc.K);
|
|
gnc.visc = Dm->Comm.sumReduce(nc.visc);
|
|
|
|
gwc.M = Dm->Comm.sumReduce(wc.M);
|
|
gwc.Px = Dm->Comm.sumReduce(wc.Px);
|
|
gwc.Py = Dm->Comm.sumReduce(wc.Py);
|
|
gwc.Pz = Dm->Comm.sumReduce(wc.Pz);
|
|
gwc.K = Dm->Comm.sumReduce(wc.K);
|
|
gwc.visc = Dm->Comm.sumReduce(wc.visc);
|
|
|
|
giwn.Mn = Dm->Comm.sumReduce(iwn.Mn);
|
|
giwn.Pnx = Dm->Comm.sumReduce(iwn.Pnx);
|
|
giwn.Pny = Dm->Comm.sumReduce(iwn.Pny);
|
|
giwn.Pnz = Dm->Comm.sumReduce(iwn.Pnz);
|
|
giwn.Kn = Dm->Comm.sumReduce(iwn.Kn);
|
|
giwn.Mw = Dm->Comm.sumReduce(iwn.Mw);
|
|
giwn.Pwx = Dm->Comm.sumReduce(iwn.Pwx);
|
|
giwn.Pwy = Dm->Comm.sumReduce(iwn.Pwy);
|
|
giwn.Pwz = Dm->Comm.sumReduce(iwn.Pwz);
|
|
giwn.Kw = Dm->Comm.sumReduce(iwn.Kw);
|
|
|
|
gifs.Mn = Dm->Comm.sumReduce(ifs.Mn);
|
|
gifs.Pnx = Dm->Comm.sumReduce(ifs.Pnx);
|
|
gifs.Pny = Dm->Comm.sumReduce(ifs.Pny);
|
|
gifs.Pnz = Dm->Comm.sumReduce(ifs.Pnz);
|
|
gifs.Mw = Dm->Comm.sumReduce(ifs.Mw);
|
|
gifs.Pwx = Dm->Comm.sumReduce(ifs.Pwx);
|
|
gifs.Pwy = Dm->Comm.sumReduce(ifs.Pwy);
|
|
gifs.Pwz = Dm->Comm.sumReduce(ifs.Pwz);
|
|
|
|
// pressure averaging
|
|
gnc.p = Dm->Comm.sumReduce(nc.p);
|
|
gnd.p = Dm->Comm.sumReduce(nd.p);
|
|
gwc.p = Dm->Comm.sumReduce(wc.p);
|
|
gwd.p = Dm->Comm.sumReduce(wd.p);
|
|
|
|
if (vol_wc_bulk > 0.0)
|
|
wc.p = wc.p / vol_wc_bulk;
|
|
if (vol_nc_bulk > 0.0)
|
|
nc.p = nc.p / vol_nc_bulk;
|
|
if (vol_wd_bulk > 0.0)
|
|
wd.p = wd.p / vol_wd_bulk;
|
|
if (vol_nd_bulk > 0.0)
|
|
nd.p = nd.p / vol_nd_bulk;
|
|
|
|
vol_wc_bulk = Dm->Comm.sumReduce(vol_wc_bulk);
|
|
vol_wd_bulk = Dm->Comm.sumReduce(vol_wd_bulk);
|
|
vol_nc_bulk = Dm->Comm.sumReduce(vol_nc_bulk);
|
|
vol_nd_bulk = Dm->Comm.sumReduce(vol_nd_bulk);
|
|
|
|
if (vol_wc_bulk > 0.0)
|
|
gwc.p = gwc.p / vol_wc_bulk;
|
|
if (vol_nc_bulk > 0.0)
|
|
gnc.p = gnc.p / vol_nc_bulk;
|
|
if (vol_wd_bulk > 0.0)
|
|
gwd.p = gwd.p / vol_wd_bulk;
|
|
if (vol_nd_bulk > 0.0)
|
|
gnd.p = gnd.p / vol_nd_bulk;
|
|
}
|
|
|
|
void SubPhase::AggregateLabels(const std::string &filename) {
|
|
|
|
int nx = Dm->Nx;
|
|
int ny = Dm->Ny;
|
|
int nz = Dm->Nz;
|
|
|
|
// assign the ID from the phase indicator field
|
|
for (int k = 0; k < nz; k++) {
|
|
for (int j = 0; j < ny; j++) {
|
|
for (int i = 0; i < nx; i++) {
|
|
int n = k * nx * ny + j * nx + i;
|
|
signed char local_id_val = Dm->id[n];
|
|
if (local_id_val > 0) {
|
|
double value = Phi(i, j, k);
|
|
if (value > 0.0)
|
|
local_id_val = 1;
|
|
else
|
|
local_id_val = 2;
|
|
}
|
|
Dm->id[n] = local_id_val;
|
|
}
|
|
}
|
|
}
|
|
Dm->Comm.barrier();
|
|
|
|
Dm->AggregateLabels(filename);
|
|
}
|