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

added analysis tools for TwoPhase object

Browse Source
This commit is contained in:
James McClure
2021-12-13 16:04:20 -05:00
parent 93f99057bb
commit 011b2f8a87
5 changed files with 528 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

288
analysis/TwoPhase.cpp
View File

@@ -14,22 +14,6 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
/*
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/TwoPhase.h"
#include "analysis/pmmc.h"
@@ -41,6 +25,8 @@
#include "IO/MeshDatabase.h"
#include "IO/Reader.h"
#include "IO/Writer.h"
#include "analysis/filters.h"
#include <memory>
@@ -430,28 +416,38 @@ void TwoPhase::UpdateSolid() {
void TwoPhase::UpdateMeshValues() {
int i, j, k, n;
fillHalo<double> fillData(Dm->Comm, Dm->rank_info, {Nx-2,Ny-2,Nz-2}, {1, 1, 1}, 0, 1);
//...........................................................................
Dm->CommunicateMeshHalo(SDn);
//Dm->CommunicateMeshHalo(SDn);
fillData.fill(SDn);
//...........................................................................
// Compute the gradients of the phase indicator and signed distance fields
pmmc_MeshGradient(SDn, SDn_x, SDn_y, SDn_z, Nx, Ny, Nz);
//...........................................................................
// Gradient of the phase indicator field
fillData.fill(SDn_x);
fillData.fill(SDn_y);
fillData.fill(SDn_z);
fillData.fill(SDs);
//...........................................................................
Dm->CommunicateMeshHalo(SDn_x);
//Dm->CommunicateMeshHalo(SDn_x);
//...........................................................................
Dm->CommunicateMeshHalo(SDn_y);
//Dm->CommunicateMeshHalo(SDn_y);
//...........................................................................
Dm->CommunicateMeshHalo(SDn_z);
//Dm->CommunicateMeshHalo(SDn_z);
//...........................................................................
Dm->CommunicateMeshHalo(SDs);
//Dm->CommunicateMeshHalo(SDs);
pmmc_MeshGradient(SDs, SDs_x, SDs_y, SDs_z, Nx, Ny, Nz);
//...........................................................................
Dm->CommunicateMeshHalo(SDs_x);
fillData.fill(SDs_x);
fillData.fill(SDs_y);
fillData.fill(SDs_z);
//Dm->CommunicateMeshHalo(SDs_x);
//...........................................................................
Dm->CommunicateMeshHalo(SDs_y);
//Dm->CommunicateMeshHalo(SDs_y);
//...........................................................................
Dm->CommunicateMeshHalo(SDs_z);
//Dm->CommunicateMeshHalo(SDs_z);
//...........................................................................
// Compute the mesh curvature of the phase indicator field
pmmc_MeshCurvature(SDn, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
@@ -579,6 +575,7 @@ void TwoPhase::ComputeLocal() {
Kwn += pmmc_CubeSurfaceInterpValue(
CubeValues, GaussCurvature, nw_pts, nw_tris, Values, i,
j, k, n_nw_pts, n_nw_tris);
Jwn += pmmc_CubeSurfaceInterpValue(
CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i,
j, k, n_nw_pts, n_nw_tris);
@@ -609,7 +606,7 @@ void TwoPhase::ComputeLocal() {
efawns += pmmc_CubeContactAngle(
CubeValues, Values, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
wwnsdnwn += pmmc_CommonCurveSpeed(
CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z, SDs_x,
SDs_y, SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
@@ -715,6 +712,192 @@ void TwoPhase::ComputeLocal() {
nonwet_morph->ComputeScalar(phase_distance, 0.f);
//printf("rank=%i completed \n",Dm->rank());
}
void TwoPhase::ComputeStatic() {
int i, j, k, n, imin, jmin, kmin, kmax;
int cube[8][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0},
{0, 0, 1}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}};
kmin = 1;
kmax = Nz - 1;
imin = jmin = 1;
/* set fluid isovalue to "grow" NWP for contact angle measurement */
fluid_isovalue = -1.5;
FILE *ANGLES = fopen("ContactAngle.csv", "a");
fprintf(ANGLES,"x y z angle\n");
for (k = kmin; k < kmax; k++) {
for (j = jmin; j < Ny - 1; j++) {
for (i = imin; i < Nx - 1; i++) {
//...........................................................................
n_nw_pts = n_ns_pts = n_ws_pts = n_nws_pts = n_local_sol_pts =
n_local_nws_pts = 0;
n_nw_tris = n_ns_tris = n_ws_tris = n_nws_seg =
n_local_sol_tris = 0;
//...........................................................................
// Compute volume averages
for (int p = 0; p < 8; p++) {
n = i + cube[p][0] + (j + cube[p][1]) * Nx +
(k + cube[p][2]) * Nx * Ny;
if (Dm->id[n] > 0) {
// 1-D index for this cube corner
// compute the norm of the gradient of the phase indicator field
// Compute the non-wetting phase volume contribution
if (Phase(i + cube[p][0], j + cube[p][1],
k + cube[p][2]) > 0) {
nwp_volume += 0.125;
} else {
wp_volume += 0.125;
}
}
}
//...........................................................................
// Construct the interfaces and common curve
pmmc_ConstructLocalCube(
SDs, SDn, solid_isovalue, fluid_isovalue, nw_pts, nw_tris,
Values, ns_pts, ns_tris, ws_pts, ws_tris, local_nws_pts,
nws_pts, nws_seg, local_sol_pts, local_sol_tris,
n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts,
n_nws_pts, n_nws_seg, i, j, k, Nx, Ny, Nz);
// wn interface averages
if (n_nw_pts > 0) {
awn += pmmc_CubeSurfaceOrientation(Gwn, nw_pts, nw_tris,
n_nw_tris);
Kwn += pmmc_CubeSurfaceInterpValue(
CubeValues, GaussCurvature, nw_pts, nw_tris, Values, i,
j, k, n_nw_pts, n_nw_tris);
Jwn += pmmc_CubeSurfaceInterpValue(
CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i,
j, k, n_nw_pts, n_nw_tris);
// Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
pmmc_CubeTrimSurfaceInterpValues(
CubeValues, MeanCurvature, SDs, nw_pts, nw_tris, Values,
DistanceValues, i, j, k, n_nw_pts, n_nw_tris, trimdist,
dummy, trJwn);
pmmc_CubeTrimSurfaceInterpInverseValues(
CubeValues, MeanCurvature, SDs, nw_pts, nw_tris, Values,
DistanceValues, i, j, k, n_nw_pts, n_nw_tris, trimdist,
dummy, trRwn);
}
// wns common curve averages
if (n_local_nws_pts > 0) {
efawns += pmmc_CubeContactAngle(
CubeValues, Values, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
for (int p = 0; p < n_local_nws_pts; p++) {
// Extract the line segment
Point A = local_nws_pts(p);
double value = Values(p);
fprintf(ANGLES, "%.8g %.8g %.8g %.8g\n", A.x, A.y, A.z, value);
}
pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x,
SDs_y, SDs_z, KNwns_values,
KGwns_values, KNwns, KGwns, nws_pts,
n_nws_pts, i, j, k);
lwns +=
pmmc_CubeCurveLength(local_nws_pts, n_local_nws_pts);
}
// Solid interface averagees
if (n_local_sol_tris > 0) {
As += pmmc_CubeSurfaceArea(local_sol_pts, local_sol_tris,
n_local_sol_tris);
// Compute the surface orientation and the interfacial area
ans += pmmc_CubeSurfaceOrientation(Gns, ns_pts, ns_tris,
n_ns_tris);
aws += pmmc_CubeSurfaceOrientation(Gws, ws_pts, ws_tris,
n_ws_tris);
}
//...........................................................................
// Compute the integral curvature of the non-wetting phase
n_nw_pts = n_nw_tris = 0;
// Compute the non-wetting phase surface and associated area
An +=
geomavg_MarchingCubes(SDn, fluid_isovalue, i, j, k, nw_pts,
n_nw_pts, nw_tris, n_nw_tris);
// Compute the integral of mean curvature
if (n_nw_pts > 0) {
pmmc_CubeTrimSurfaceInterpValues(
CubeValues, MeanCurvature, SDs, nw_pts, nw_tris, Values,
DistanceValues, i, j, k, n_nw_pts, n_nw_tris, trimdist,
trawn, dummy);
}
Jn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature,
nw_pts, nw_tris, Values, i, j,
k, n_nw_pts, n_nw_tris);
// Compute Euler characteristic from integral of gaussian curvature
Kn += pmmc_CubeSurfaceInterpValue(CubeValues, GaussCurvature,
nw_pts, nw_tris, Values, i, j,
k, n_nw_pts, n_nw_tris);
euler += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
n_nw_tris, i, j, k);
}
}
}
fclose(ANGLES);
Array<char> phase_label(Nx, Ny, Nz);
Array<double> phase_distance(Nx, Ny, Nz);
// Analyze the wetting fluid
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 (!(Dm->id[n] > 0)) {
// Solid phase
phase_label(i, j, k) = 1;
} else if (SDn(i, j, k) < 0.0) {
// wetting phase
phase_label(i, j, k) = 0;
} else {
// non-wetting phase
phase_label(i, j, k) = 1;
}
phase_distance(i, j, k) =
2.0 * double(phase_label(i, j, k)) - 1.0;
}
}
}
CalcDist(phase_distance, phase_label, *Dm);
wet_morph->ComputeScalar(phase_distance, 0.f);
//printf("generating distance at rank=%i \n",Dm->rank());
// Analyze the wetting fluid
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 (!(Dm->id[n] > 0)) {
// Solid phase
phase_label(i, j, k) = 1;
} else if (SDn(i, j, k) < 0.0) {
// wetting phase
phase_label(i, j, k) = 1;
} else {
// non-wetting phase
phase_label(i, j, k) = 0;
}
phase_distance(i, j, k) =
2.0 * double(phase_label(i, j, k)) - 1.0;
}
}
}
CalcDist(phase_distance, phase_label, *Dm);
nonwet_morph->ComputeScalar(phase_distance, 0.f);
}
void TwoPhase::AssignComponentLabels() {
//int LabelNWP=1;
@@ -1283,6 +1466,7 @@ void TwoPhase::Reduce() {
van_global(2) = van_global(2) / nwp_volume_global;
}
// Normalize surface averages by the interfacial area
/*
if (awn_global > 0.0) {
Jwn_global /= awn_global;
Kwn_global /= awn_global;
@@ -1299,6 +1483,7 @@ void TwoPhase::Reduce() {
for (i = 0; i < 3; i++)
vawns_global(i) /= lwns_global;
}
*/
if (trawn_global > 0.0) {
trJwn_global /= trawn_global;
trRwn_global /= trawn_global;
@@ -1314,15 +1499,17 @@ void TwoPhase::Reduce() {
Gws_global(i) /= aws_global;
euler_global /= (2 * PI);
//sat_w = 1.0 - nwp_volume_global*iVol_global/porosity;
sat_w = 1.0 - nwp_volume_global / (nwp_volume_global + wp_volume_global);
// Compute the specific interfacial areas and common line length (dimensionless per unit volume)
/*
awn_global = awn_global * iVol_global;
ans_global = ans_global * iVol_global;
aws_global = aws_global * iVol_global;
dEs = dEs * iVol_global;
lwns_global = lwns_global * iVol_global;
*/
}
void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) {
@@ -1334,6 +1521,53 @@ void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) {
lwns_global *= D * D;
}
void TwoPhase::PrintStatic() {
if (Dm->rank() == 0) {
FILE *STATIC;
STATIC = fopen("geometry.csv", "a+");
if (fseek(STATIC, 0, SEEK_SET) == fseek(STATIC, 0, SEEK_CUR)) {
// If timelog is empty, write a short header to list the averages
fprintf(STATIC, "sw awn ans aws Jwn Kwn lwns cwns KGws "
"KGwn "); // Scalar averages
fprintf(STATIC,
"Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
fprintf(STATIC, "Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
fprintf(STATIC, "Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
fprintf(STATIC, "trawn trJwn trRwn "); //trimmed curvature,
fprintf(STATIC, "Vw Aw Jw Xw "); //miknowski measures,
fprintf(STATIC, "Vn An Jn Xn\n"); //miknowski measures,
//fprintf(STATIC,"Euler Kn2 Jn2 An2\n"); //miknowski measures,
}
fprintf(STATIC, "%.5g ", sat_w); // saturation
fprintf(STATIC, "%.5g %.5g %.5g ", awn_global, ans_global,
aws_global); // interfacial areas
fprintf(STATIC, "%.5g %.5g ", Jwn_global,
Kwn_global); // curvature of wn interface
fprintf(STATIC, "%.5g ", lwns_global); // common curve length
fprintf(STATIC, "%.5g ", efawns_global); // average contact angle
fprintf(STATIC, "%.5g %.5g ", KNwns_global,
KGwns_global); // total curvature contributions of common line
fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gwn_global(0),
Gwn_global(1), Gwn_global(2), Gwn_global(3), Gwn_global(4),
Gwn_global(5)); // orientation of wn interface
fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gns_global(0),
Gns_global(1), Gns_global(2), Gns_global(3), Gns_global(4),
Gns_global(5)); // orientation of ns interface
fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gws_global(0),
Gws_global(1), Gws_global(2), Gws_global(3), Gws_global(4),
Gws_global(5)); // orientation of ws interface
fprintf(STATIC, "%.5g %.5g %.5g ", trawn_global, trJwn_global,
trRwn_global); // Trimmed curvature
fprintf(STATIC, "%.5g %.5g %.5g %.5g ", wet_morph->V(), wet_morph->A(),
wet_morph->H(), wet_morph->X());
fprintf(STATIC, "%.5g %.5g %.5g %.5g\n", nonwet_morph->V(),
nonwet_morph->A(), nonwet_morph->H(), nonwet_morph->X());
//fprintf(STATIC,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global); // minkowski measures
fclose(STATIC);
}
}
void TwoPhase::PrintAll(int timestep) {
if (Dm->rank() == 0) {
fprintf(TIMELOG, "%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",

2
analysis/TwoPhase.h
View File

@@ -184,6 +184,8 @@ public:
void ColorToSignedDistance(double Beta, DoubleArray &ColorData,
DoubleArray &DistData);
void ComputeLocal();
void ComputeStatic();
void PrintStatic();
void AssignComponentLabels();
void ComponentAverages();
void Reduce();

74
analysis/pmmc.h
View File

@@ -4057,7 +4057,7 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues,
(A.z - B.z) * (A.z - B.z));
integral += 0.5 * length * (vA + vB);
}
return integral;
}
//--------------------------------------------------------------------------------------------------------
@@ -4420,7 +4420,9 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
double twnsx, twnsy, twnsz, nwnsx, nwnsy, nwnsz,
K; // tangent,normal and curvature
double nsx, nsy, nsz, norm;
double nwx, nwy, nwz;
double nwsx, nwsy, nwsz;
double nwnx, nwny, nwnz;
Point P, A, B;
// Local trilinear approximation for tangent and normal vector
@@ -4430,47 +4432,18 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
for (k = kc; k < kc + 2; k++) {
for (j = jc; j < jc + 2; j++) {
for (i = ic; i < ic + 2; i++) {
// Compute all of the derivatives using finite differences
// fluid phase indicator field
// fx = 0.5*(f(i+1,j,k) - f(i-1,j,k));
// fy = 0.5*(f(i,j+1,k) - f(i,j-1,k));
// fz = 0.5*(f(i,j,k+1) - f(i,j,k-1));
/*fxx = f(i+1,j,k) - 2.0*f(i,j,k) + f(i-1,j,k);
fyy = f(i,j+1,k) - 2.0*f(i,j,k) + f(i,j-1,k);
fzz = f(i,j,k+1) - 2.0*f(i,j,k) + f(i,j,k-1);
fxy = 0.25*(f(i+1,j+1,k) - f(i+1,j-1,k) - f(i-1,j+1,k) + f(i-1,j-1,k));
fxz = 0.25*(f(i+1,j,k+1) - f(i+1,j,k-1) - f(i-1,j,k+1) + f(i-1,j,k-1));
fyz = 0.25*(f(i,j+1,k+1) - f(i,j+1,k-1) - f(i,j-1,k+1) + f(i,j-1,k-1));
*/
// solid distance function
// sx = 0.5*(s(i+1,j,k) - s(i-1,j,k));
// sy = 0.5*(s(i,j+1,k) - s(i,j-1,k));
// sz = 0.5*(s(i,j,k+1) - s(i,j,k-1));
/* sxx = s(i+1,j,k) - 2.0*s(i,j,k) + s(i-1,j,k);
syy = s(i,j+1,k) - 2.0*s(i,j,k) + s(i,j-1,k);
szz = s(i,j,k+1) - 2.0*s(i,j,k) + s(i,j,k-1);
sxy = 0.25*(s(i+1,j+1,k) - s(i+1,j-1,k) - s(i-1,j+1,k) + s(i-1,j-1,k));
sxz = 0.25*(s(i+1,j,k+1) - s(i+1,j,k-1) - s(i-1,j,k+1) + s(i-1,j,k-1));
syz = 0.25*(s(i,j+1,k+1) - s(i,j+1,k-1) - s(i,j-1,k+1) + s(i,j-1,k-1));
*/
/* // Compute the Jacobean matrix for tangent vector
Axx = sxy*fz + sy*fxz - sxz*fy - sz*fxy;
Axy = sxz*fx + sz*fxx - sxx*fz - sx*fxz;
Axz = sxx*fy + sx*fxy - sxy*fx - sy*fxx;
Ayx = syy*fz + sy*fyz - syz*fy - sz*fyy;
Ayy = syz*fx + sz*fxy - sxy*fz - sx*fyz;
Ayz = sxy*fy + sx*fyy - syy*fx - sy*fxy;
Azx = syz*fz + sy*fzz - szz*fy - sz*fyz;
Azy = szz*fx + sz*fxz - sxz*fz - sx*fzz;
Azz = sxz*fy + sx*fyz - syz*fx - sy*fxz;
*/
sx = s_x(i, j, k);
sy = s_y(i, j, k);
sz = s_z(i, j, k);
fx = f_x(i, j, k);
fy = f_y(i, j, k);
fz = f_z(i, j, k);
// Normal to fluid surface
Nx.Corners(i - ic, j - jc, k - kc) = fx;
Ny.Corners(i - ic, j - jc, k - kc) = fy;
Nz.Corners(i - ic, j - jc, k - kc) = fz;
// Normal to solid surface
Sx.Corners(i - ic, j - jc, k - kc) = sx;
Sy.Corners(i - ic, j - jc, k - kc) = sy;
@@ -4577,8 +4550,19 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
nsx /= norm;
nsy /= norm;
nsz /= norm;
// Normal vector to the fluid surface
nwx = Nx.eval(P);
nwy = Ny.eval(P);
nwz = Nz.eval(P);
norm = sqrt(nwx * nwx + nwy * nwy + nwz * nwz);
if (norm == 0.0)
norm = 1.0;
nwx /= norm;
nwy /= norm;
nwz /= norm;
// normal in the surface tangent plane (rel. geodesic curvature)
// common curve normal in the solid surface tangent plane (rel. geodesic curvature)
nwsx = twnsy * nsz - twnsz * nsy;
nwsy = twnsz * nsx - twnsx * nsz;
nwsz = twnsx * nsy - twnsy * nsx;
@@ -4588,16 +4572,28 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
nwsx /= norm;
nwsy /= norm;
nwsz /= norm;
if (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz < 0.0) {
nwnsx = -nwnsx;
nwnsy = -nwnsy;
nwnsz = -nwnsz;
}
// common curve normal in the fluid surface tangent plane (rel. geodesic curvature)
nwnx = twnsy * nwz - twnsz * nwy;
nwny = twnsz * nwx - twnsx * nwz;
nwnz = twnsx * nwy - twnsy * nwx;
norm = sqrt(nwnx * nwnx + nwny * nwny + nwnz * nwnz);
if (norm == 0.0)
norm = 1.0;
nwnx /= norm;
nwny /= norm;
nwnz /= norm;
if (length > 0.0) {
// normal curvature component in the direction of the solid surface
KNavg += K * (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz) * length;
//KNavg += K * (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz) * length;
KNavg += K * (nwnx * nwnsx + nwny * nwnsy + nwnz * nwnsz) * length;
//geodesic curvature
KGavg += K * (nwsx * nwnsx + nwsy * nwnsy + nwsz * nwnsz) * length;
}

1
tests/CMakeLists.txt
View File

@@ -39,6 +39,7 @@ ADD_LBPM_EXECUTABLE( convertIO )
ADD_LBPM_EXECUTABLE( DataAggregator )
#ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel )(
ADD_LBPM_EXECUTABLE( lbpm_minkowski_scalar )
ADD_LBPM_EXECUTABLE( lbpm_TwoPhase_analysis )
ADD_LBPM_EXECUTABLE( TestPoissonSolver )
ADD_LBPM_EXECUTABLE( TestIonModel )
ADD_LBPM_EXECUTABLE( TestNernstPlanck )

229
tests/lbpm_TwoPhase_analysis.cpp Normal file
View File

@@ -0,0 +1,229 @@
/*
* Compute porous media marching cubes analysis (PMMC) based on input image data
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <functional>
#include "common/Array.h"
#include "common/Domain.h"
#include "common/Communication.h"
#include "common/MPI.h"
#include "IO/MeshDatabase.h"
#include "IO/Mesh.h"
#include "IO/Writer.h"
#include "IO/netcdf.h"
#include "analysis/analysis.h"
#include "analysis/filters.h"
#include "analysis/distance.h"
#include "analysis/Minkowski.h"
#include "analysis/TwoPhase.h"
#include "ProfilerApp.h"
int main(int argc, char **argv)
{
// Initialize MPI
Utilities::startup( argc, argv );
Utilities::MPI comm( MPI_COMM_WORLD );
int rank = comm.getRank();
//int nprocs = comm.getSize();
{
Utilities::setErrorHandlers();
PROFILE_START("Main");
//std::vector<std::string> filenames;
if ( argc<2 ) {
if ( rank == 0 ){
printf("At least one filename must be specified\n");
}
return 1;
}
std::string filename = std::string(argv[1]);
if ( rank == 0 ){
printf("Input data file: %s\n",filename.c_str());
}
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
auto vis_db = db->getDatabase( "Visualization" );
// Read domain parameters
auto Filename = domain_db->getScalar<std::string>( "Filename" );
auto size = domain_db->getVector<int>( "n" );
auto SIZE = domain_db->getVector<int>( "N" );
auto nproc = domain_db->getVector<int>( "nproc" );
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
auto nx = size[0];
auto ny = size[1];
auto nz = size[2];
int i,j,k,n;
std::shared_ptr<Domain> Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis
comm.barrier();
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;
// Initialize the object
Dm->id[n] = 1;
}
}
}
Dm->CommInit();
/* read the data */
if (domain_db->keyExists( "Filename" )){
auto Filename = domain_db->getScalar<std::string>( "Filename" );
Dm->Decomp(Filename);
}
else{
Dm->ReadIDs();
}
fillHalo<double> fillDouble(Dm->Comm, Dm->rank_info, {nx,ny,nz}, {1, 1, 1}, 0, 1);
// Compute the Minkowski functionals
comm.barrier();
std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
std::shared_ptr<Minkowski> Geometry(new Minkowski(Dm));
// Calculate the distance
// Initialize the domain and communication
nx+=2; ny+=2; nz+=2;
Array<char> id(nx,ny,nz);
DoubleArray Distance(nx,ny,nz);
/* * * * * * * * * * * * * * * * * * * * * */
// Solve for the position of the solid phase
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;
// Initialize the object
if (Dm->id[n] == ReadValues[0]) id(i,j,k) = 0;
else id(i,j,k) = 1;
}
}
}
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;
// Initialize distance to +/- 1
Averages->SDs(i,j,k) = 2.0*double(id(i,j,k))-1.0;
}
}
}
fillDouble.fill(Averages->SDs);
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(Averages->SDs,id,*Dm);
/* * * * * * * * * * * * * * * * * * * * * */
// Solve for the position of the fluid phase
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;
// Initialize the object
if (Dm->id[n] == ReadValues[1]){
id(i,j,k) = 1;
Averages->Phase(i,j,k) = 1.0;
}
else {
id(i,j,k) = 0;
Averages->Phase(i,j,k) = -1.0;
}
}
}
}
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;
// Initialize distance to +/- 1
Distance(i,j,k) = 2.0*double(id(i,j,k))-1.0;
}
}
}
fillDouble.fill(Averages->Phase);
fillDouble.fill(Distance);
if (rank==0) printf("Initialized fluid phase -- Converting to Signed Distance function \n");
CalcDist(Distance,id,*Dm);
Mean3D(Distance,Averages->SDn);
/* * * * * * * * * * * * * * * * * * * * * */
if (rank==0) printf("Computing Minkowski functionals \n");
Geometry->ComputeScalar(Distance,0.f);
Geometry->PrintAll();
Averages->Initialize();
Averages->UpdateMeshValues();
Averages->ComputeStatic();
Averages->Reduce();
Averages->PrintStatic();
/*
Averages.Initialize();
Averages.ComponentAverages();
Averages.SortBlobs();
Averages.PrintComponents(timestep);
*/
auto format = vis_db->getWithDefault<string>("format", "hdf5");
vis_db = db->getDatabase("Visualization");
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);
auto PhaseVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
IO::initialize("", format, "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);
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VariableType::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2);
visData[0].vars.push_back(SignDistVar);
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);
Array<double> &SignData = visData[0].vars[0]->data;
Array<double> &PhaseData = visData[0].vars[1]->data;
ASSERT(visData[0].vars[0]->name == "SignDist");
ASSERT(visData[0].vars[1]->name == "Phase");
fillData.copy(Averages->SDs, SignData);
fillData.copy(Averages->SDn, PhaseData);
int timestep = 0;
IO::writeData(timestep, visData, Dm->Comm);
}
PROFILE_STOP("Main");
PROFILE_SAVE("TwoPhase",true);
Utilities::shutdown();
return 0;
}