fix whitespace merge issue
This commit is contained in:
James E McClure
commit
1f6d37208e
|
|
@ -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>
|
||||
|
||||
|
|
@ -375,6 +361,7 @@ void TwoPhase::Initialize() {
|
|||
wwndnw = 0.0;
|
||||
wwnsdnwn = 0.0;
|
||||
Jwnwwndnw = 0.0;
|
||||
Xwn = Xns = Xws = 0.0;
|
||||
}
|
||||
|
||||
void TwoPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB,
|
||||
|
|
@ -430,28 +417,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 +576,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 +607,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 +713,217 @@ 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.0;
|
||||
|
||||
string FILENAME = "ContactAngle";
|
||||
|
||||
char LocalRankString[8];
|
||||
char LocalRankFilename[40];
|
||||
sprintf(LocalRankString, "%05d", Dm->rank());
|
||||
sprintf(LocalRankFilename, "%s%s%s", "ContactAngle.", LocalRankString,".csv");
|
||||
FILE *ANGLES = fopen(LocalRankFilename, "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);
|
||||
|
||||
Xwn += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
|
||||
n_nw_tris, i, j, k);
|
||||
|
||||
// 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);
|
||||
|
||||
/* half contribution for vertices / edges at the common line
|
||||
* each cube with contact line has a net of undercounting vertices
|
||||
* each cube is undercounting edges due to internal counts
|
||||
*/
|
||||
Xwn += 0.25*n_local_nws_pts - 0.5;
|
||||
Xws += 0.25*n_local_nws_pts - 0.5;
|
||||
Xns += 0.25*n_local_nws_pts - 0.5;
|
||||
}
|
||||
|
||||
// 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);
|
||||
|
||||
Xws += geomavg_EulerCharacteristic(ws_pts, ws_tris, n_ws_pts,
|
||||
n_ws_tris, i, j, k);
|
||||
|
||||
Xns += geomavg_EulerCharacteristic(ns_pts, ns_tris, n_ns_pts,
|
||||
n_ns_tris, i, j, k);
|
||||
}
|
||||
//...........................................................................
|
||||
// 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;
|
||||
|
|
@ -1219,7 +1428,7 @@ void TwoPhase::ComponentAverages() {
|
|||
|
||||
void TwoPhase::Reduce() {
|
||||
int i;
|
||||
double iVol_global = 1.0 / Volume;
|
||||
//double iVol_global = 1.0 / Volume;
|
||||
//...........................................................................
|
||||
Dm->Comm.barrier();
|
||||
nwp_volume_global = Dm->Comm.sumReduce(nwp_volume);
|
||||
|
|
@ -1258,10 +1467,14 @@ void TwoPhase::Reduce() {
|
|||
trawn_global = Dm->Comm.sumReduce(trawn);
|
||||
trJwn_global = Dm->Comm.sumReduce(trJwn);
|
||||
trRwn_global = Dm->Comm.sumReduce(trRwn);
|
||||
euler_global = Dm->Comm.sumReduce(euler);
|
||||
Xwn_global = Dm->Comm.sumReduce(Xwn);
|
||||
Xws_global = Dm->Comm.sumReduce(Xws);
|
||||
Xns_global = Dm->Comm.sumReduce(Xns);
|
||||
An_global = Dm->Comm.sumReduce(An);
|
||||
Jn_global = Dm->Comm.sumReduce(Jn);
|
||||
Kn_global = Dm->Comm.sumReduce(Kn);
|
||||
euler_global = Dm->Comm.sumReduce(euler);
|
||||
|
||||
Dm->Comm.barrier();
|
||||
|
||||
// Normalize the phase averages
|
||||
|
|
@ -1285,17 +1498,18 @@ void TwoPhase::Reduce() {
|
|||
// Normalize surface averages by the interfacial area
|
||||
if (awn_global > 0.0) {
|
||||
Jwn_global /= awn_global;
|
||||
Kwn_global /= awn_global;
|
||||
//Kwn_global /= awn_global;
|
||||
wwndnw_global /= awn_global;
|
||||
for (i = 0; i < 3; i++)
|
||||
vawn_global(i) /= awn_global;
|
||||
for (i = 0; i < 6; i++)
|
||||
Gwn_global(i) /= awn_global;
|
||||
}
|
||||
|
||||
if (lwns_global > 0.0) {
|
||||
efawns_global /= lwns_global;
|
||||
KNwns_global /= lwns_global;
|
||||
KGwns_global /= lwns_global;
|
||||
//KNwns_global /= lwns_global;
|
||||
//KGwns_global /= lwns_global;
|
||||
for (i = 0; i < 3; i++)
|
||||
vawns_global(i) /= lwns_global;
|
||||
}
|
||||
|
|
@ -1314,15 +1528,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 +1550,55 @@ 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 Xwn Xws Xns "); // 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 ", Xwn_global, Xns_global,
|
||||
Xws_global); // Euler characteristic
|
||||
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 ",
|
||||
|
|
|
|||
|
|
@ -103,7 +103,9 @@ public:
|
|||
double lwns_global;
|
||||
double efawns, efawns_global; // averaged contact angle
|
||||
double euler, Kn, Jn, An;
|
||||
double Xwn, Xns, Xws;
|
||||
double euler_global, Kn_global, Jn_global, An_global;
|
||||
double Xwn_global, Xns_global, Xws_global;
|
||||
|
||||
double rho_n, rho_w;
|
||||
double nu_n, nu_w;
|
||||
|
|
@ -184,6 +186,8 @@ public:
|
|||
void ColorToSignedDistance(double Beta, DoubleArray &ColorData,
|
||||
DoubleArray &DistData);
|
||||
void ComputeLocal();
|
||||
void ComputeStatic();
|
||||
void PrintStatic();
|
||||
void AssignComponentLabels();
|
||||
void ComponentAverages();
|
||||
void Reduce();
|
||||
|
|
|
|||
|
|
@ -137,8 +137,7 @@ void Morphology::Initialize(std::shared_ptr<Domain> Dm, DoubleArray &Distance) {
|
|||
morphRadius.resize(recvLoc);
|
||||
//..............................
|
||||
/* send the morphological radius */
|
||||
Dm->Comm.Irecv(&morphRadius[recvOffset_x], recvCount, Dm->rank_x(),
|
||||
recvtag + 0);
|
||||
Dm->Comm.Irecv(&morphRadius[recvOffset_x], recvCount, Dm->rank_x(), recvtag + 0);
|
||||
Dm->Comm.send(&tmpDistance[0], sendCount, Dm->rank_X(), sendtag + 0);
|
||||
/* send the shift values */
|
||||
Dm->Comm.Irecv(&xShift[recvOffset_x], recvCount, Dm->rank_x(), recvtag + 1);
|
||||
|
|
@ -502,7 +501,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
|
|||
if (rank == 0)
|
||||
printf("Maximum pore size: %f \n", maxdistGlobal);
|
||||
final_void_fraction = volume_fraction; //initialize
|
||||
|
||||
|
||||
int ii, jj, kk;
|
||||
int imin, jmin, kmin, imax, jmax, kmax;
|
||||
int Nx = nx;
|
||||
|
|
@ -524,27 +523,28 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
|
|||
|
||||
int numTry = 0;
|
||||
int maxTry = 100;
|
||||
while (void_fraction_new > VoidFraction && numTry < maxTry) {
|
||||
while ( !(void_fraction_new < VoidFraction) && numTry < maxTry) {
|
||||
numTry++;
|
||||
void_fraction_diff_old = void_fraction_diff_new;
|
||||
void_fraction_old = void_fraction_new;
|
||||
Rcrit_old = Rcrit_new;
|
||||
Rcrit_new -= deltaR * Rcrit_old;
|
||||
if (rank==0) printf("Try %i with radius %f \n", numTry, Rcrit_new);
|
||||
if (Rcrit_new < 0.5) {
|
||||
numTry = maxTry;
|
||||
}
|
||||
int Window = round(Rcrit_new);
|
||||
if (Window == 0)
|
||||
Window =
|
||||
1; // If Window = 0 at the begining, after the following process will have sw=1.0
|
||||
Window = 1; // If Window = 0 at the begining, after the following process will have sw=1.0
|
||||
// and sw<Sw will be immediately broken
|
||||
double LocalNumber = 0.f;
|
||||
for (int k = 0; k < Nz; k++) {
|
||||
for (int j = 0; j < Ny; j++) {
|
||||
for (int i = 0; i < Nx; i++) {
|
||||
for (int k = 1; k < Nz-1; k++) {
|
||||
for (int j = 1; j < Ny-1; j++) {
|
||||
for (int i = 1; i < Nx-1; i++) {
|
||||
n = k * nx * ny + j * nx + i;
|
||||
if (SignDist(i, j, k) > Rcrit_new) {
|
||||
// loop over the window and update
|
||||
//printf("Distance(%i %i %i) = %f \n",i,j,k, SignDist(i,j,k));
|
||||
imin = max(1, i - Window);
|
||||
jmin = max(1, j - Window);
|
||||
kmin = max(1, k - Window);
|
||||
|
|
@ -571,7 +571,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
|
|||
}
|
||||
}
|
||||
}
|
||||
LocalNumber += Structure.GetOverlaps(Dm, id, ErodeLabel, NewLabel);
|
||||
//LocalNumber += Structure.GetOverlaps(Dm, id, ErodeLabel, NewLabel);
|
||||
|
||||
count = 0.f;
|
||||
for (int k = 1; k < Nz - 1; k++) {
|
||||
|
|
@ -611,7 +611,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
|
|||
|
||||
//***************************************************************************************
|
||||
double MorphDrain(DoubleArray &SignDist, signed char *id,
|
||||
std::shared_ptr<Domain> Dm, double VoidFraction) {
|
||||
std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius) {
|
||||
// SignDist is the distance to the object that you want to constaing the morphological opening
|
||||
// VoidFraction is the the empty space where the object inst
|
||||
// id is a labeled map
|
||||
|
|
@ -688,6 +688,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id,
|
|||
double deltaR = 0.05; // amount to change the radius in voxel units
|
||||
double Rcrit_old = maxdistGlobal;
|
||||
double Rcrit_new = maxdistGlobal;
|
||||
|
||||
if (InitialRadius < maxdistGlobal){
|
||||
Rcrit_old = InitialRadius;
|
||||
Rcrit_new = InitialRadius;
|
||||
}
|
||||
//if (argc>2){
|
||||
// Rcrit_new = strtod(argv[2],NULL);
|
||||
// if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new);
|
||||
|
|
|
|||
|
|
@ -7,7 +7,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
|
|||
std::shared_ptr<Domain> Dm, double VoidFraction,
|
||||
signed char ErodeLabel, signed char ReplaceLabel);
|
||||
double MorphDrain(DoubleArray &SignDist, signed char *id,
|
||||
std::shared_ptr<Domain> Dm, double VoidFraction);
|
||||
std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius);
|
||||
double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
|
||||
std::shared_ptr<Domain> Dm, double TargetVol,
|
||||
double WallFactor);
|
||||
|
|
|
|||
|
|
@ -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,35 @@ 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;
|
||||
/* normal to ws interface boundary should point into fluid (same direction as gradient) */
|
||||
if (nwx * nwsx + nwy * nwsy + nwz * nwsz < 0.0) {
|
||||
nwsx = -nwsx;
|
||||
nwsy = -nwsy;
|
||||
nwsz = -nwsz;
|
||||
}
|
||||
|
||||
// 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;
|
||||
/* normal to wn interface boundary should point into the solid */
|
||||
if (nsx * nwnx + nsy * nwny + nsz * nwnz > 0.0) {
|
||||
nwnx = -nwnx;
|
||||
nwny = -nwny;
|
||||
nwnz = -nwnz;
|
||||
}
|
||||
|
||||
|
||||
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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -296,7 +296,7 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity
|
|||
* @param finish - lattice node to finish loop
|
||||
* @param Np - size of local sub-domain (derived from Domain structure)
|
||||
*/
|
||||
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,
|
||||
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB, bool UseSlippingVelBC,
|
||||
int start, int finish, int Np);
|
||||
|
||||
/**
|
||||
|
|
@ -307,19 +307,22 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *di
|
|||
* @param Psi -
|
||||
* @param ElectricField - electric field
|
||||
* @param tau - relaxation time
|
||||
* @param epsilon_LB -
|
||||
* @param epsilon_LB - dielectric constant of medium
|
||||
* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve
|
||||
* @param start - lattice node to start loop
|
||||
* @param finish - lattice node to finish loop
|
||||
* @param Np - size of local sub-domain (derived from Domain structure)
|
||||
*/
|
||||
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau,
|
||||
double epsilon_LB, int start, int finish, int Np);
|
||||
double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np);
|
||||
/**
|
||||
* \brief Poisson-Boltzmann solver / solve electric potential based on AA odd access pattern for D3Q7
|
||||
* @param neighborList - neighbors based on D3Q19 lattice structure
|
||||
* @param Map - mapping between sparse and dense representations
|
||||
* @param dist - D3Q7 distributions
|
||||
* @param Psi -
|
||||
* @param epsilon_LB - dielectric constant of medium
|
||||
* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve
|
||||
* @param start - lattice node to start loop
|
||||
* @param finish - lattice node to finish loop
|
||||
* @param Np - size of local sub-domain (derived from Domain structure)
|
||||
|
|
@ -350,10 +353,10 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in
|
|||
|
||||
// LBM Stokes Model (adapted from MRT model)
|
||||
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
|
||||
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np);
|
||||
double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np);
|
||||
|
||||
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB,
|
||||
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np);
|
||||
double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np);
|
||||
|
||||
extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np);
|
||||
|
||||
|
|
|
|||
|
|
@ -95,7 +95,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(
|
|||
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
|
||||
double *dist, double *Den_charge,
|
||||
double *Psi, double *ElectricField,
|
||||
double tau, double epsilon_LB,
|
||||
double tau, double epsilon_LB, bool UseSlippingVelBC,
|
||||
int start, int finish, int Np) {
|
||||
int n;
|
||||
double psi; //electric potential
|
||||
|
|
@ -109,8 +109,9 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
|
|||
for (n = start; n < finish; n++) {
|
||||
|
||||
//Load data
|
||||
rho_e = Den_charge[n];
|
||||
rho_e = rho_e / epsilon_LB;
|
||||
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
|
||||
//and thus the net space charge density is zero.
|
||||
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
|
||||
idx = Map[n];
|
||||
psi = Psi[idx];
|
||||
|
||||
|
|
@ -175,8 +176,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
|
|||
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
|
||||
double *Den_charge, double *Psi,
|
||||
double *ElectricField, double tau,
|
||||
double epsilon_LB, int start,
|
||||
int finish, int Np) {
|
||||
double epsilon_LB, bool UseSlippingVelBC,
|
||||
int start, int finish, int Np) {
|
||||
int n;
|
||||
double psi; //electric potential
|
||||
double Ex, Ey, Ez; //electric field
|
||||
|
|
@ -188,8 +189,9 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
|
|||
for (n = start; n < finish; n++) {
|
||||
|
||||
//Load data
|
||||
rho_e = Den_charge[n];
|
||||
rho_e = rho_e / epsilon_LB;
|
||||
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
|
||||
//and thus the net space charge density is zero.
|
||||
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
|
||||
idx = Map[n];
|
||||
psi = Psi[idx];
|
||||
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
|
|||
double *dist, double *Velocity, double *ChargeDensity,
|
||||
double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
|
||||
double Gy, double Gz, double rho0, double den_scale, double h,
|
||||
double time_conv, int start, int finish, int Np) {
|
||||
double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
|
||||
double fq;
|
||||
// conserved momemnts
|
||||
double rho, jx, jy, jz;
|
||||
|
|
@ -38,13 +38,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
|
|||
Ey = ElectricField[n + 1 * Np];
|
||||
Ez = ElectricField[n + 2 * Np];
|
||||
//compute total body force, including input body force (Gx,Gy,Gz)
|
||||
Fx =
|
||||
Gx +
|
||||
rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
|
||||
// q=0
|
||||
|
|
@ -479,7 +477,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
|
|||
int *neighborList, double *dist, double *Velocity, double *ChargeDensity,
|
||||
double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
|
||||
double Gy, double Gz, double rho0, double den_scale, double h,
|
||||
double time_conv, int start, int finish, int Np) {
|
||||
double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
|
||||
double fq;
|
||||
// conserved momemnts
|
||||
double rho, jx, jy, jz;
|
||||
|
|
@ -513,12 +511,21 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
|
|||
Ex = ElectricField[n + 0 * Np];
|
||||
Ey = ElectricField[n + 1 * Np];
|
||||
Ez = ElectricField[n + 2 * Np];
|
||||
|
||||
//compute total body force, including input body force (Gx,Gy,Gz)
|
||||
Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
//Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
// den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
//Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
// den_scale;
|
||||
//Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
// den_scale;
|
||||
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
|
||||
//and body force induced by external efectric field is reduced to slipping velocity BC.
|
||||
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
|
||||
// q=0
|
||||
|
|
|
|||
|
|
@ -104,7 +104,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub
|
|||
}
|
||||
}
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double psi;//electric potential
|
||||
|
|
@ -122,8 +122,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
|
|||
if (n<finish) {
|
||||
|
||||
//Load data
|
||||
rho_e = Den_charge[n];
|
||||
rho_e = rho_e/epsilon_LB;
|
||||
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
|
||||
//and thus the net space charge density is zero.
|
||||
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
|
||||
idx=Map[n];
|
||||
psi = Psi[idx];
|
||||
|
||||
|
|
@ -184,7 +185,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
|
|||
}
|
||||
}
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double psi;//electric potential
|
||||
|
|
@ -201,8 +202,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *
|
|||
if (n<finish) {
|
||||
|
||||
//Load data
|
||||
rho_e = Den_charge[n];
|
||||
rho_e = rho_e/epsilon_LB;
|
||||
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
|
||||
//and thus the net space charge density is zero.
|
||||
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
|
||||
idx=Map[n];
|
||||
psi = Psi[idx];
|
||||
|
||||
|
|
@ -293,10 +295,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
|
|||
//cudaProfilerStop();
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
|
||||
|
||||
cudaError_t err = cudaGetLastError();
|
||||
if (cudaSuccess != err){
|
||||
|
|
@ -305,10 +307,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
|
|||
//cudaProfilerStop();
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
|
||||
|
||||
cudaError_t err = cudaGetLastError();
|
||||
if (cudaSuccess != err){
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@
|
|||
#define NBLOCKS 1024
|
||||
#define NTHREADS 256
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double fq;
|
||||
|
|
@ -46,9 +46,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
|
|||
Ey = ElectricField[n+1*Np];
|
||||
Ez = ElectricField[n+2*Np];
|
||||
//compute total body force, including input body force (Gx,Gy,Gz)
|
||||
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
|
||||
// q=0
|
||||
fq = dist[n];
|
||||
|
|
@ -510,7 +513,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
|
|||
}
|
||||
}
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double fq;
|
||||
|
|
@ -550,9 +553,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
|
|||
Ey = ElectricField[n+1*Np];
|
||||
Ez = ElectricField[n+2*Np];
|
||||
//compute total body force, including input body force (Gx,Gy,Gz)
|
||||
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
|
||||
// q=0
|
||||
fq = dist[n];
|
||||
|
|
@ -969,10 +975,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
|
|||
}
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np);
|
||||
|
||||
cudaError_t err = cudaGetLastError();
|
||||
if (cudaSuccess != err){
|
||||
|
|
@ -981,10 +987,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
|
|||
//cudaProfilerStop();
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np);
|
||||
|
||||
cudaError_t err = cudaGetLastError();
|
||||
if (cudaSuccess != err){
|
||||
|
|
|
|||
137
docs/source/examples/analysis/TwoPhase.rst
Normal file
137
docs/source/examples/analysis/TwoPhase.rst
Normal file
|
|
@ -0,0 +1,137 @@
|
|||
*************************
|
||||
Measuring Contact Angles
|
||||
*************************
|
||||
|
||||
LBPM includes specialized data analysis capabilities for two-fluid systems. While these
|
||||
components are generally designed for in situ analysis of simulation data, they can also
|
||||
be applied independently to analyze 3D image data. In this example we consider applying
|
||||
the analysis tools implemented in ``lbpm_TwoPhase_analysis``, which are designed to
|
||||
analyze two-fluid configurations in porous media. The numerical implementation used to
|
||||
construct the common line are described in ( https://doi.org/10.1016/j.advwatres.2006.06.010 ).
|
||||
Methods used to measure the contact angle are described in ( https://doi.org/10.1017/jfm.2016.212 ).
|
||||
|
||||
Source files for the example are included in the LBPM repository
|
||||
in the directory ``examples/Droplet``. A simple python code is included
|
||||
to set up a fluid droplet on a flat surface
|
||||
|
||||
.. code:: python
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pylab as plt
|
||||
|
||||
D=np.ones((80,80,40),dtype="uint8")
|
||||
|
||||
cx = 40
|
||||
cy = 40
|
||||
cz = 0
|
||||
|
||||
for i in range(0,80):
|
||||
for j in range (0, 80):
|
||||
for k in range (0,40):
|
||||
dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
|
||||
if (dist < 25) :
|
||||
D[i,j,k] = 2
|
||||
if k < 4 :
|
||||
D[i,j,k] = 0
|
||||
|
||||
D.tofile("droplet_40x80x80.raw")
|
||||
|
||||
|
||||
The input file provided below will specify how the analysis should be performed. The name and the dimensions of the input
|
||||
file are provided in the ``Domain`` section, as with other LBPM simulators. For large images, additional processors can be
|
||||
used to speed up the analysis or take advantage of distributed memory. The ``ReadValues`` list should be used to specify the
|
||||
labels to use for analysis. The first label will be taken to be the solid. The second label will be taken to be the fluid
|
||||
to analyze, which in this case will be the droplet labeled with ``2`` above.
|
||||
|
||||
.. code:: bash
|
||||
|
||||
Domain {
|
||||
Filename = "droplet_40x80x80.raw"
|
||||
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
|
||||
n = 40, 80, 80 // Size of local domain (Nx,Ny,Nz)
|
||||
N = 40, 80, 80 // size of the input image
|
||||
voxel_length = 1.0
|
||||
BC = 0 // Boundary condition type
|
||||
ReadType = "8bit"
|
||||
ReadValues = 0, 2, 1
|
||||
WriteValues = 0, 2, 1
|
||||
}
|
||||
Visualization {
|
||||
}
|
||||
|
||||
|
||||
The analysis can be launched as ``mpirun -np 1 $LBPM_DIR/lbpm_TwoPhase_analysis input.db``. Output should appear as
|
||||
follows:
|
||||
|
||||
.. code:: bash
|
||||
|
||||
Input data file: input.db
|
||||
voxel length = 1.000000 micron
|
||||
Input media: droplet_40x80x80.raw
|
||||
Relabeling 3 values
|
||||
oldvalue=0, newvalue =0
|
||||
oldvalue=2, newvalue =2
|
||||
oldvalue=1, newvalue =1
|
||||
Dimensions of segmented image: 40 x 80 x 80
|
||||
Reading 8-bit input data
|
||||
Read segmented data from droplet_40x80x80.raw
|
||||
Label=0, Count=25600
|
||||
Label=2, Count=25773
|
||||
Label=1, Count=204627
|
||||
Distributing subdomains across 1 processors
|
||||
Process grid: 1 x 1 x 1
|
||||
Subdomain size: 40 x 80 x 80
|
||||
Size of transition region: 0
|
||||
Media porosity = 0.900000
|
||||
Initialized solid phase -- Converting to Signed Distance function
|
||||
Initialized fluid phase -- Converting to Signed Distance function
|
||||
Computing Minkowski functionals
|
||||
|
||||
The ``TwoPhase`` analysis class will generate signed distance functions for the solid and fluid surfaces.
|
||||
Using the distance functions, the interfaces and common line will be constructed. Contact angles are logged
|
||||
for each processor, e.g. ``ContactAngle.00000.csv``, which specifies the x, y, z coordinates for each measurement
|
||||
along with the cosine of the contact angle. Averaged measures (determined over the entire input image)
|
||||
are logged to ``geometry.csv``
|
||||
|
||||
* ``sw`` -- water saturation
|
||||
* ``awn`` -- surface area of meniscus between wn fluids
|
||||
* ``ans`` -- surface area between fluid n and solid
|
||||
* ``aws`` -- surface area between fluid w and solid
|
||||
* ``Jwn`` -- integral of mean curvature of meniscus
|
||||
* ``Kwn`` -- integral of Gaussian curvature of meniscus
|
||||
* ``lwns`` -- length of common line
|
||||
* ``cwns`` -- average contact angle
|
||||
* ``KGws`` -- geodesic curvature of common line relative to ws surface
|
||||
* ``KGwn`` -- geodesic curvature of common line relative to wn surface
|
||||
* ``Gwnxx`` -- orientation tensor component for wn surface
|
||||
* ``Gwnyy`` -- orientation tensor component for wn surface
|
||||
* ``Gwnzz`` -- orientation tensor component for wn surface
|
||||
* ``Gwnxy`` -- orientation tensor component for wn surface
|
||||
* ``Gwnxz`` -- orientation tensor component for wn surface
|
||||
* ``Gwnyz`` -- orientation tensor component for wn surface
|
||||
* ``Gwsxx`` -- orientation tensor component for ws surface
|
||||
* ``Gwsyy`` -- orientation tensor component for ws surface
|
||||
* ``Gwszz`` -- orientation tensor component for ws surface
|
||||
* ``Gwsxy`` -- orientation tensor component for ws surface
|
||||
* ``Gwsxz`` -- orientation tensor component for ws surface
|
||||
* ``Gwsyz`` -- orientation tensor component for ws surface
|
||||
* ``Gnsxx`` -- orientation tensor component for ns surface
|
||||
* ``Gnsyy`` -- orientation tensor component for ns surface
|
||||
* ``Gnszz`` -- orientation tensor component for ns surface
|
||||
* ``Gnsxy`` -- orientation tensor component for ns surface
|
||||
* ``Gnsxz`` -- orientation tensor component for ns surface
|
||||
* ``Gnsyz`` -- orientation tensor component for ns surface
|
||||
* ``trawn`` -- trimmed surface area for meniscus (one voxel from solid)
|
||||
* ``trJwn`` -- mean curvature for trimmed meniscus
|
||||
* ``trRwn`` -- radius of curvature for trimmed meniscus
|
||||
* ``Vw`` -- volume of fluid w
|
||||
* ``Aw`` -- boundary surface area for fluid w
|
||||
* ``Jw`` -- integral of mean curvature for fluid w
|
||||
* ``Xw`` -- Euler characteristic for fluid w
|
||||
* ``Vn`` -- volume of fluid n
|
||||
* ``An`` -- boundary surface area for fluid n
|
||||
* ``Jn`` -- integral of mean curvature for fluid n
|
||||
* ``Xn`` -- Euler characteristic for fluid n
|
||||
|
||||
|
||||
|
||||
125
docs/source/examples/color/fractionalFlow.rst
Normal file
125
docs/source/examples/color/fractionalFlow.rst
Normal file
|
|
@ -0,0 +1,125 @@
|
|||
********************************
|
||||
Steady-state fractional flow
|
||||
********************************
|
||||
|
||||
In this example we simulate a steady-state flow with a constant driving force. This will enforce a periodic boundary condition
|
||||
in all directions. While the driving force may be set in any direction, we will set it in the z-direction to be consistent
|
||||
with the convention for pressure and velocity boundary conditions.
|
||||
|
||||
|
||||
For the case considered in ``example/DiscPack`` we specify the following information in the input file
|
||||
|
||||
.. code:: c
|
||||
|
||||
Domain {
|
||||
Filename = "discs_3x128x128.raw.morphdrain.raw"
|
||||
ReadType = "8bit" // data type
|
||||
N = 3, 128, 128 // size of original image
|
||||
nproc = 1, 2, 2 // process grid
|
||||
n = 3, 64, 64 // sub-domain size
|
||||
voxel_length = 1.0 // voxel length (in microns)
|
||||
ReadValues = 0, 1, 2 // labels within the original image
|
||||
WriteValues = 0, 1, 2 // associated labels to be used by LBPM
|
||||
BC = 0 // fully periodic BC
|
||||
Sw = 0.35 // target saturation for morphological tools
|
||||
}
|
||||
|
||||
Color {
|
||||
protocol = "fractional flow"
|
||||
capillary_number = -1e-5 // capillary number for the displacement, positive="oil injection"
|
||||
timestepMax = 200000 // maximum timtestep
|
||||
alpha = 0.01 // controls interfacial tension
|
||||
rhoA = 1.0 // controls the density of fluid A
|
||||
rhoB = 1.0 // controls the density of fluid B
|
||||
tauA = 0.7 // controls the viscosity of fluid A
|
||||
tauB = 0.7 // controls the viscosity of fluid B
|
||||
F = 0, 0, 1e-5 // body force
|
||||
WettingConvention = "SCAL"
|
||||
ComponentLabels = 0 // image labels for solid voxels
|
||||
ComponentAffinity = 0.9 // controls the wetting affinity for each label
|
||||
Restart = false
|
||||
}
|
||||
Analysis {
|
||||
analysis_interval = 1000 // logging interval for timelog.csv
|
||||
subphase_analysis_interval = 500000 // loggging interval for subphase.csv
|
||||
N_threads = 4 // number of analysis threads (GPU version only)
|
||||
visualization_interval = 10000 // interval to write visualization files
|
||||
restart_interval = 10000000 // interval to write restart file
|
||||
restart_file = "Restart" // base name of restart file
|
||||
}
|
||||
Visualization {
|
||||
format = "hdf5"
|
||||
write_silo = true // write SILO databases with assigned variables
|
||||
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
|
||||
save_phase_field = true // save phase field within SILO database
|
||||
save_pressure = true // save pressure field within SILO database
|
||||
save_velocity = false // save velocity field within SILO database
|
||||
}
|
||||
FlowAdaptor {
|
||||
max_steady_timesteps = 25000 // maximum number of timesteps per steady point
|
||||
min_steady_timesteps = 25000 // minimum number of timesteps per steady point
|
||||
fractional_flow_increment = 0.0003 // parameter that controls rate of mass seeding
|
||||
skip_timesteps = 10000 // number of timesteps to spend in flow adaptor
|
||||
}
|
||||
|
||||
|
||||
Once this has been set, we launch ``lbpm_color_simulator`` in the same way as other parallel tools
|
||||
|
||||
.. code:: bash
|
||||
|
||||
mpirun -np 4 $LBPM_BIN/lbpm_color_simulator input.db
|
||||
|
||||
Successful output looks like the following
|
||||
|
||||
|
||||
.. code:: bash
|
||||
|
||||
********************************************************
|
||||
Running Color LBM
|
||||
********************************************************
|
||||
voxel length = 1.000000 micron
|
||||
voxel length = 1.000000 micron
|
||||
Input media: discs_3x128x128.raw.morphdrain.raw
|
||||
Relabeling 3 values
|
||||
oldvalue=0, newvalue =0
|
||||
oldvalue=1, newvalue =1
|
||||
oldvalue=2, newvalue =2
|
||||
Dimensions of segmented image: 3 x 128 x 128
|
||||
Reading 8-bit input data
|
||||
Read segmented data from discs_3x128x128.raw.morphdrain.raw
|
||||
Label=0, Count=11862
|
||||
Label=1, Count=26430
|
||||
Label=2, Count=10860
|
||||
Distributing subdomains across 4 processors
|
||||
Process grid: 1 x 2 x 2
|
||||
Subdomain size: 3 x 64 x 64
|
||||
Size of transition region: 0
|
||||
Media porosity = 0.758667
|
||||
Initialized solid phase -- Converting to Signed Distance function
|
||||
Domain set.
|
||||
Create ScaLBL_Communicator
|
||||
Set up memory efficient layout, 9090 | 9120 | 21780
|
||||
Allocating distributions
|
||||
Setting up device map and neighbor list
|
||||
Component labels: 1
|
||||
label=0, affinity=-0.900000, volume fraction==0.417582
|
||||
Initializing distributions
|
||||
Initializing phase field
|
||||
Affinities - rank 0:
|
||||
Main: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 2: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 3: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 4: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Affinities - rank 0:
|
||||
Main: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 2: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 3: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
Thread 4: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
|
||||
********************************************************
|
||||
CPU time = 0.001501
|
||||
Lattice update rate (per core)= 6.074861 MLUPS
|
||||
Lattice update rate (per MPI process)= 6.074861 MLUPS
|
||||
(flatten density field)
|
||||
|
||||
|
|
@ -11,51 +11,57 @@ For the case considered in ``example/DiscPack`` we specify the following informa
|
|||
|
||||
.. code:: c
|
||||
|
||||
Domain {
|
||||
Filename = "discs_3x128x128.raw.morphdrain.raw"
|
||||
ReadType = "8bit" // data type
|
||||
N = 3, 128, 128 // size of original image
|
||||
nproc = 1, 2, 2 // process grid
|
||||
n = 3, 64, 64 // sub-domain size
|
||||
voxel_length = 1.0 // voxel length (in microns)
|
||||
ReadValues = 0, 1, 2 // labels within the original image
|
||||
WriteValues = 0, 1, 2 // associated labels to be used by LBPM
|
||||
BC = 0 // fully periodic BC
|
||||
Sw = 0.35 // target saturation for morphological tools
|
||||
}
|
||||
Domain {
|
||||
Filename = "discs_3x128x128.raw.morphdrain.raw"
|
||||
ReadType = "8bit" // data type
|
||||
N = 3, 128, 128 // size of original image
|
||||
nproc = 1, 2, 2 // process grid
|
||||
n = 3, 64, 64 // sub-domain size
|
||||
voxel_length = 1.0 // voxel length (in microns)
|
||||
ReadValues = 0, 1, 2 // labels within the original image
|
||||
WriteValues = 0, 1, 2 // associated labels to be used by LBPM
|
||||
BC = 0 // fully periodic BC
|
||||
Sw = 0.35 // target saturation for morphological tools
|
||||
}
|
||||
|
||||
Color {
|
||||
capillary_number = -1e-5 // capillary number for the displacement, positive="oil injection"
|
||||
timestepMax = 20000 // maximum timtestep
|
||||
alpha = 0.01 // controls interfacial tension
|
||||
rhoA = 1.0 // controls the density of fluid A
|
||||
rhoB = 1.0 // controls the density of fluid B
|
||||
tauA = 0.7 // controls the viscosity of fluid A
|
||||
tauB = 0.7 // controls the viscosity of fluid B
|
||||
F = 0, 0, 1e-5 // body force
|
||||
WettingConvention = "SCAL"
|
||||
ComponentLabels = 0 // image labels for solid voxels
|
||||
ComponentAffinity = 0.9 // controls the wetting affinity for each label
|
||||
Restart = false
|
||||
}
|
||||
Analysis {
|
||||
analysis_interval = 1000 // logging interval for timelog.csv
|
||||
subphase_analysis_interval = 500000 // loggging interval for subphase.csv
|
||||
N_threads = 4 // number of analysis threads (GPU version only)
|
||||
visualization_interval = 10000 // interval to write visualization files
|
||||
restart_interval = 10000000 // interval to write restart file
|
||||
restart_file = "Restart" // base name of restart file
|
||||
}
|
||||
Visualization {
|
||||
format = "hdf5"
|
||||
write_silo = true // write SILO databases with assigned variables
|
||||
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
|
||||
save_phase_field = true // save phase field within SILO database
|
||||
save_pressure = true // save pressure field within SILO database
|
||||
save_velocity = false // save velocity field within SILO database
|
||||
}
|
||||
FlowAdaptor {
|
||||
}
|
||||
Color {
|
||||
protocol = "fractional flow"
|
||||
capillary_number = 1e-5 // capillary number for the displacement, positive="oil injection"
|
||||
timestepMax = 500000 // maximum timtestep
|
||||
alpha = 0.005 // controls interfacial tension
|
||||
rhoA = 1.0 // controls the density of fluid A
|
||||
rhoB = 1.0 // controls the density of fluid B
|
||||
tauA = 0.7 // controls the viscosity of fluid A
|
||||
tauB = 0.7 // controls the viscosity of fluid B
|
||||
F = 0, 0, 1e-5 // body force
|
||||
WettingConvention = "SCAL"
|
||||
ComponentLabels = 0 // image labels for solid voxels
|
||||
ComponentAffinity = 0.9 // controls the wetting affinity for each label
|
||||
Restart = false
|
||||
}
|
||||
Analysis {
|
||||
analysis_interval = 1000 // logging interval for timelog.csv
|
||||
subphase_analysis_interval = 500000 // loggging interval for subphase.csv
|
||||
N_threads = 4 // number of analysis threads (GPU version only)
|
||||
visualization_interval = 1000000 // interval to write visualization files
|
||||
restart_interval = 10000000 // interval to write restart file
|
||||
restart_file = "Restart" // base name of restart file
|
||||
}
|
||||
Visualization {
|
||||
format = "hdf5"
|
||||
write_silo = true // write SILO databases with assigned variables
|
||||
save_8bit_raw = true // write labeled 8-bit binary files with phase assignments
|
||||
save_phase_field = true // save phase field within SILO database
|
||||
save_pressure = true // save pressure field within SILO database
|
||||
save_velocity = false // save velocity field within SILO database
|
||||
}
|
||||
FlowAdaptor {
|
||||
min_steady_timesteps = 250000 // minimum number of timesteps per steady point
|
||||
max_steady_timesteps = 300000 // maximum number of timesteps per steady point
|
||||
fractional_flow_increment = 0.1 // parameter that controls rate of mass seeding
|
||||
skip_timesteps = 10000 // number of timesteps to spend in flow adaptor
|
||||
endpoint_threshold = 0.1 // endpoint exit criterion
|
||||
}
|
||||
|
||||
|
||||
Once this has been set, we launch ``lbpm_color_simulator`` in the same way as other parallel tools
|
||||
|
|
|
|||
|
|
@ -18,4 +18,5 @@ a basic introduction to working with LBPM.
|
|||
morphology/*
|
||||
|
||||
color/*
|
||||
|
||||
|
||||
analysis/*
|
||||
|
|
|
|||
|
|
@ -1,13 +1,15 @@
|
|||
============
|
||||
\============
|
||||
Publications
|
||||
============
|
||||
|
||||
* James E McClure, Zhe Li, Mark Berrill, Thomas Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9
|
||||
* J.E. McClure, Z. Li, M. Berrill, T. Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9
|
||||
|
||||
|
||||
* James E. McClure, Zhe Li, Adrian P. Sheppard, Cass T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670
|
||||
* J. E. McClure, Z. Li, A.P. Sheppard, C.T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670
|
||||
|
||||
|
||||
* Y.D. Wang, T. Chung, R.T. Armstrong, J. McClure, T. Ramstad, P. Mostaghimi, "Accelerated Computation of Relative Permeability by Coupled Morphological and Direct Multiphase Flow Simulation" Journal of Computational Physics (401) (2020) https://doi.org/10.1016/j.jcp.2019.108966
|
||||
|
||||
|
||||
* J.E. McClure, M. Berrill, W. Gray, C.T. Miller, C.T. "Tracking interface and common curve dynamics for two-fluid flow in porous media. Journal of Fluid Mechanics, 796, 211-232 (2016) https://doi.org/10.1017/jfm.2016.212
|
||||
|
||||
* J.E.McClure, D.Adalsteinsson, C.Pan, W.G.Gray, C.T.Miller "Approximation of interfacial properties in multiphase porous medium systems" Advances in Water Resources, 30 (3): 354-365 (2007) https://doi.org/10.1016/j.advwatres.2006.06.010
|
||||
|
|
|
|||
|
|
@ -120,7 +120,7 @@ two fluids are permitted to freely mix between the endpoints. Beyond the endpoin
|
|||
term is used to drive spontaneous imbibition into the grey voxels
|
||||
|
||||
|
||||
..math::
|
||||
.. math::
|
||||
:nowrap:
|
||||
|
||||
$$
|
||||
|
|
|
|||
|
|
@ -3,6 +3,8 @@
|
|||
INSTALL_EXAMPLE( Bubble )
|
||||
INSTALL_EXAMPLE( ConstrainedBubble )
|
||||
INSTALL_EXAMPLE( Piston )
|
||||
INSTALL_EXAMPLE( Droplet )
|
||||
INSTALL_EXAMPLE( DropletCoalescence )
|
||||
INSTALL_EXAMPLE( Plates )
|
||||
INSTALL_EXAMPLE( SquareTube )
|
||||
INSTALL_EXAMPLE( InkBottle )
|
||||
|
|
|
|||
197
example/Droplet/CreateDroplet.ipynb
Normal file
197
example/Droplet/CreateDroplet.ipynb
Normal file
File diff suppressed because one or more lines are too long
20
example/Droplet/CreateDroplet.py
Normal file
20
example/Droplet/CreateDroplet.py
Normal file
|
|
@ -0,0 +1,20 @@
|
|||
import numpy as np
|
||||
import matplotlib.pylab as plt
|
||||
|
||||
D=np.ones((80,80,40),dtype="uint8")
|
||||
|
||||
cx = 40
|
||||
cy = 40
|
||||
cz = 0
|
||||
|
||||
for i in range(0,80):
|
||||
for j in range (0, 80):
|
||||
for k in range (0,40):
|
||||
dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
|
||||
if (dist < 25) :
|
||||
D[i,j,k] = 2
|
||||
if k < 4 :
|
||||
D[i,j,k] = 0
|
||||
|
||||
D.tofile("droplet_40x80x80.raw")
|
||||
|
||||
14
example/Droplet/input.db
Normal file
14
example/Droplet/input.db
Normal file
|
|
@ -0,0 +1,14 @@
|
|||
Domain {
|
||||
Filename = "droplet_40x80x80.raw"
|
||||
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
|
||||
n = 40, 80, 80 // Size of local domain (Nx,Ny,Nz)
|
||||
N = 40, 80, 80 // size of the input image
|
||||
voxel_length = 1.0
|
||||
BC = 0 // Boundary condition type
|
||||
ReadType = "8bit"
|
||||
ReadValues = 0, 2, 1
|
||||
WriteValues = 0, 2, 1
|
||||
}
|
||||
|
||||
Visualization {
|
||||
}
|
||||
|
|
@ -104,7 +104,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub
|
|||
}
|
||||
}
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double psi;//electric potential
|
||||
|
|
@ -122,8 +122,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
|
|||
if (n<finish) {
|
||||
|
||||
//Load data
|
||||
rho_e = Den_charge[n];
|
||||
rho_e = rho_e/epsilon_LB;
|
||||
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
|
||||
//and thus the net space charge density is zero.
|
||||
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
|
||||
idx=Map[n];
|
||||
psi = Psi[idx];
|
||||
|
||||
|
|
@ -184,7 +185,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
|
|||
}
|
||||
}
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double psi;//electric potential
|
||||
|
|
@ -201,8 +202,9 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *
|
|||
if (n<finish) {
|
||||
|
||||
//Load data
|
||||
rho_e = Den_charge[n];
|
||||
rho_e = rho_e/epsilon_LB;
|
||||
//When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
|
||||
//and thus the net space charge density is zero.
|
||||
rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
|
||||
idx=Map[n];
|
||||
psi = Psi[idx];
|
||||
|
||||
|
|
@ -293,10 +295,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
|
|||
//cudaProfilerStop();
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
|
||||
|
||||
hipError_t err = hipGetLastError();
|
||||
if (hipSuccess != err){
|
||||
|
|
@ -305,10 +307,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
|
|||
//cudaProfilerStop();
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
|
||||
|
||||
hipError_t err = hipGetLastError();
|
||||
if (hipSuccess != err){
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@
|
|||
#define NTHREADS 256
|
||||
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double fq;
|
||||
|
|
@ -47,9 +47,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
|
|||
Ey = ElectricField[n+1*Np];
|
||||
Ez = ElectricField[n+2*Np];
|
||||
//compute total body force, including input body force (Gx,Gy,Gz)
|
||||
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
|
||||
// q=0
|
||||
fq = dist[n];
|
||||
|
|
@ -511,7 +514,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
|
|||
}
|
||||
}
|
||||
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
|
||||
__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC,int start, int finish, int Np){
|
||||
|
||||
int n;
|
||||
double fq;
|
||||
|
|
@ -551,9 +554,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
|
|||
Ey = ElectricField[n+1*Np];
|
||||
Ez = ElectricField[n+2*Np];
|
||||
//compute total body force, including input body force (Gx,Gy,Gz)
|
||||
Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
|
||||
Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale; //the extra factors at the end necessarily convert unit from phys to LB
|
||||
Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
|
||||
den_scale;
|
||||
|
||||
// q=0
|
||||
fq = dist[n];
|
||||
|
|
@ -970,10 +976,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
|
|||
}
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
|
||||
|
||||
hipError_t err = hipGetLastError();
|
||||
if (hipSuccess != err){
|
||||
|
|
@ -982,10 +988,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
|
|||
//cudaProfilerStop();
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
|
||||
extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
|
||||
|
||||
//cudaProfilerStart();
|
||||
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
|
||||
dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
|
||||
|
||||
hipError_t err = hipGetLastError();
|
||||
if (hipSuccess != err){
|
||||
|
|
|
|||
|
|
@ -885,6 +885,14 @@ double ScaLBL_ColorModel::Run(int returntime) {
|
|||
double vBd_x = Averages->gwd.Px / Mass_w;
|
||||
double vBd_y = Averages->gwd.Py / Mass_w;
|
||||
double vBd_z = Averages->gwd.Pz / Mass_w;
|
||||
|
||||
/* contribution from water films */
|
||||
double water_film_flow_rate =
|
||||
Averages->gwb.V * (Averages->giwn.Pwx * dir_x + Averages->giwn.Pwy * dir_y + Averages->giwn.Pwz * dir_z) /
|
||||
Averages->gwb.M / Dm->Volume;
|
||||
double not_water_film_flow_rate =
|
||||
Averages->gnb.V * (Averages->giwn.Pnx * dir_x + Averages->giwn.Pny * dir_y + Averages->giwn.Pnz * dir_z) /
|
||||
Averages->gnb.M / Dm->Volume;
|
||||
|
||||
double flow_rate_A_connected =
|
||||
Vol_nc *
|
||||
|
|
@ -921,6 +929,11 @@ double ScaLBL_ColorModel::Run(int returntime) {
|
|||
double kAeff = h * h * muA * (flow_rate_A) / (force_mag);
|
||||
double kBeff = h * h * muB * (flow_rate_B) / (force_mag);
|
||||
|
||||
|
||||
/* not counting films */
|
||||
double krnf = kAeff - h * h * muA * not_water_film_flow_rate / force_mag;
|
||||
double krwf = kBeff - h * h * muB * water_film_flow_rate / force_mag;
|
||||
|
||||
// Saturation normalized effective permeability to account for decoupled phases and
|
||||
// effective porosity.
|
||||
double kAeff_low = (1.0 - current_saturation) * h * h *
|
||||
|
|
@ -945,6 +958,8 @@ double ScaLBL_ColorModel::Run(int returntime) {
|
|||
fprintf(kr_log_file, "timesteps sat.water ");
|
||||
fprintf(kr_log_file, "eff.perm.oil.upper.bound "
|
||||
"eff.perm.water.upper.bound ");
|
||||
fprintf(kr_log_file, "eff.perm.oil.film "
|
||||
"eff.perm.water.film ");
|
||||
fprintf(kr_log_file, "eff.perm.oil.lower.bound "
|
||||
"eff.perm.water.lower.bound ");
|
||||
fprintf(kr_log_file,
|
||||
|
|
@ -962,6 +977,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
|
|||
fprintf(kr_log_file, "%i %.5g ", CURRENT_TIMESTEP,
|
||||
current_saturation);
|
||||
fprintf(kr_log_file, "%.5g %.5g ", kAeff, kBeff);
|
||||
fprintf(kr_log_file, "%.5g %.5g ", krnf, krwf);
|
||||
fprintf(kr_log_file, "%.5g %.5g ", kAeff_low, kBeff_low);
|
||||
fprintf(kr_log_file, "%.5g %.5g ", kAeff_connected,
|
||||
kBeff_connected);
|
||||
|
|
|
|||
|
|
@ -523,7 +523,7 @@ void ScaLBL_Poisson::Initialize(double time_conv_from_Study){
|
|||
//}
|
||||
}
|
||||
|
||||
void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
|
||||
void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study){
|
||||
|
||||
//.......create and start timer............
|
||||
//double starttime,stoptime,cputime;
|
||||
|
|
@ -537,13 +537,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
|
|||
// *************ODD TIMESTEP*************//
|
||||
timestep++;
|
||||
SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential
|
||||
SolvePoissonAAodd(ChargeDensity);//perform collision
|
||||
SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision
|
||||
ScaLBL_Comm->Barrier(); comm.barrier();
|
||||
|
||||
// *************EVEN TIMESTEP*************//
|
||||
timestep++;
|
||||
SolveElectricPotentialAAeven(timestep_from_Study);//update electric potential
|
||||
SolvePoissonAAeven(ChargeDensity);//perform collision
|
||||
SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision
|
||||
ScaLBL_Comm->Barrier(); comm.barrier();
|
||||
//************************************************************************/
|
||||
|
||||
|
|
@ -724,9 +724,9 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){
|
|||
ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
|
||||
}
|
||||
|
||||
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){
|
||||
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
|
||||
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np);
|
||||
void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC){
|
||||
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
|
||||
ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
|
||||
//TODO: perhaps add another ScaLBL_Comm routine to update Psi values on solid boundary nodes.
|
||||
//something like:
|
||||
//ScaLBL_Comm->SolidDirichletBoundaryUpdates(Psi, Psi_BCLabel, timestep);
|
||||
|
|
@ -739,9 +739,9 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){
|
|||
//}
|
||||
}
|
||||
|
||||
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity){
|
||||
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
|
||||
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np);
|
||||
void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC){
|
||||
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
|
||||
ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np);
|
||||
ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel);
|
||||
//if (BoundaryConditionSolid==1){
|
||||
// ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi);
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@ public:
|
|||
void ReadInput();
|
||||
void Create();
|
||||
void Initialize(double time_conv_from_Study);
|
||||
void Run(double *ChargeDensity, int timestep_from_Study);
|
||||
void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study);
|
||||
void getElectricPotential(DoubleArray &ReturnValues);
|
||||
void getElectricPotential_debug(int timestep);
|
||||
void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y,
|
||||
|
|
@ -111,8 +111,8 @@ private:
|
|||
void SolveElectricPotentialAAodd(int timestep_from_Study);
|
||||
void SolveElectricPotentialAAeven(int timestep_from_Study);
|
||||
//void SolveElectricField();
|
||||
void SolvePoissonAAodd(double *ChargeDensity);
|
||||
void SolvePoissonAAeven(double *ChargeDensity);
|
||||
void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC);
|
||||
void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC);
|
||||
void getConvergenceLog(int timestep,double error);
|
||||
double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step);
|
||||
|
||||
|
|
|
|||
|
|
@ -593,7 +593,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
|
|||
ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
|
||||
ScaLBL_D3Q19_AAodd_StokesMRT(
|
||||
NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA,
|
||||
rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv,
|
||||
rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
|
||||
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
|
||||
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
|
||||
// Set boundary conditions
|
||||
|
|
@ -610,8 +610,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
|
|||
}
|
||||
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity,
|
||||
ElectricField, rlx_setA, rlx_setB, Fx, Fy,
|
||||
Fz, rho0, den_scale, h, time_conv, 0,
|
||||
ScaLBL_Comm->LastExterior(), Np);
|
||||
Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
|
||||
0, ScaLBL_Comm->LastExterior(), Np);
|
||||
|
||||
if (UseSlippingVelBC == true) {
|
||||
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
|
||||
|
|
@ -626,8 +626,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
|
|||
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
|
||||
ScaLBL_D3Q19_AAeven_StokesMRT(
|
||||
fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx,
|
||||
Fy, Fz, rho0, den_scale, h, time_conv, ScaLBL_Comm->FirstInterior(),
|
||||
ScaLBL_Comm->LastInterior(), Np);
|
||||
Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
|
||||
ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
|
||||
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
|
||||
// Set boundary conditions
|
||||
if (BoundaryCondition == 3) {
|
||||
|
|
@ -643,8 +643,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
|
|||
}
|
||||
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity,
|
||||
ElectricField, rlx_setA, rlx_setB, Fx, Fy,
|
||||
Fz, rho0, den_scale, h, time_conv, 0,
|
||||
ScaLBL_Comm->LastExterior(), Np);
|
||||
Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
|
||||
0, ScaLBL_Comm->LastExterior(), Np);
|
||||
if (UseSlippingVelBC == true) {
|
||||
ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
|
||||
fq, ZetaPotentialSolid, ElectricField, SolidGrad, epsilon_LB,
|
||||
|
|
|
|||
9
sample_scripts/Hipify.sh
Normal file
9
sample_scripts/Hipify.sh
Normal file
|
|
@ -0,0 +1,9 @@
|
|||
#!/bin/bash
|
||||
|
||||
echo "Hipify cuda file $1"
|
||||
|
||||
sed -i 's/cudaError_t/hipError_t/g' $1
|
||||
sed -i 's/cudaGetLastError/hipGetLastError/g' $1
|
||||
sed -i 's/cudaSuccess/hipSuccess/g' $1
|
||||
sed -i 's/cudaGetErrorString/hipGetErrorString/g' $1
|
||||
sed -i 's/CUDA error/hip error/g' $1
|
||||
|
|
@ -5,19 +5,16 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc \
|
|||
-D CMAKE_C_FLAGS="-O3 -fPIC" \
|
||||
-D CMAKE_CXX_FLAGS="-O3 -fPIC " \
|
||||
-D CMAKE_CXX_STANDARD=14 \
|
||||
-D MPIEXEC=mpirun \
|
||||
-D MPIEXEC=mpirun \
|
||||
-D CMAKE_BUILD_PARALLEL_LEVEL=4 \
|
||||
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
|
||||
-D CMAKE_BUILD_TYPE:STRING=Release \
|
||||
-D CUDA_FLAGS="-arch sm_35" \
|
||||
-D CUDA_HOST_COMPILER="/usr/bin/gcc" \
|
||||
-D USE_HDF5=1 \
|
||||
-D HDF5_DIRECTORY="/opt/arden/hdf5/1.8.12" \
|
||||
-D HDF5_LIB="/opt/arden/hdf5/1.8.12/lib/libhdf5.a"\
|
||||
-D USE_SILO=1 \
|
||||
-D SILO_LIB="/opt/arden/silo/4.10.2/lib/libsiloh5.a" \
|
||||
-D SILO_DIRECTORY="/opt/arden/silo/4.10.2" \
|
||||
-D USE_NETCDF=0 \
|
||||
-D NETCDF_DIRECTORY="/opt/arden/netcdf/4.6.1" \
|
||||
-D USE_CUDA=0 \
|
||||
-D USE_TIMER=0 \
|
||||
~/Programs/LBPM
|
||||
~/Programs/LBPM-WIA
|
||||
|
||||
# -D HDF5_LIB="/opt/arden/hdf5/1.8.12/lib/libhdf5.a"\
|
||||
|
|
|
|||
|
|
@ -2,19 +2,23 @@
|
|||
#module load cmake
|
||||
|
||||
#source /gpfs/gpfs_stage1/b6p315aa/setup/setup-mpi.sh
|
||||
module load cmake gcc
|
||||
module load cmake gcc/7.5.0
|
||||
module load cuda
|
||||
#/ccs/proj/csc380/mcclurej
|
||||
export HDF5_DIR=/ccs/proj/csc380/mcclurej/install/hdf5/1.8.12/
|
||||
export SILO_DIR=/ccs/proj/csc380/mcclurej/install/silo/4.10.2/
|
||||
export NETCDF_DIR=/ccs/proj/geo136/install/netcdf/4.6.1
|
||||
module load hdf5
|
||||
|
||||
#/ccs/proj/csc380/mcclurej
|
||||
#export HDF5_DIR=/ccs/proj/csc380/mcclurej/install/hdf5/1.8.12/
|
||||
#export SILO_DIR=/ccs/proj/csc380/mcclurej/install/silo/4.10.2/
|
||||
#export NETCDF_DIR=/ccs/proj/geo136/install/netcdf/4.6.1
|
||||
export HDF5_DIR="$OLCF_HDF5_ROOT"
|
||||
# configure
|
||||
|
||||
rm -rf CMake*
|
||||
cmake \
|
||||
-D CMAKE_BUILD_TYPE:STRING=Release \
|
||||
-D CMAKE_C_COMPILER:PATH=mpicc \
|
||||
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
|
||||
-D CMAKE_CXX_FLAGS="-lstdc++" \
|
||||
-D CMAKE_CXX_STANDARD=14 \
|
||||
-D USE_CUDA=1 \
|
||||
-D CMAKE_CUDA_FLAGS="-arch sm_70 -Xptxas=-v -Xptxas -dlcm=cg -lineinfo" \
|
||||
|
|
@ -25,7 +29,7 @@ cmake \
|
|||
-D USE_HDF5=1 \
|
||||
-D HDF5_DIRECTORY="$HDF5_DIR" \
|
||||
-D HDF5_LIB="$HDF5_DIR/lib/libhdf5.a" \
|
||||
-D USE_SILO=1 \
|
||||
-D USE_SILO=0 \
|
||||
-D SILO_LIB="$SILO_DIR/lib/libsiloh5.a" \
|
||||
-D SILO_DIRECTORY="$SILO_DIR" \
|
||||
-D USE_NETCDF=0 \
|
||||
|
|
|
|||
33
sample_scripts/configure_tinker_A100
Executable file
33
sample_scripts/configure_tinker_A100
Executable file
|
|
@ -0,0 +1,33 @@
|
|||
# load the module for cmake
|
||||
module reset
|
||||
module load CMake
|
||||
#module load CUDA
|
||||
#module load OpenMPI
|
||||
#module load cuda11.2/toolkit
|
||||
|
||||
module load OpenMPI/4.0.5-gcccuda-2020b
|
||||
module load HDF5
|
||||
|
||||
#export GCC_BIN=/apps/easybuild/software/tinkercliffs-rome_a100/GCCcore/10.2.0/bin
|
||||
|
||||
cmake \
|
||||
-D CMAKE_C_COMPILER:PATH=mpicc \
|
||||
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
|
||||
-D CMAKE_C_FLAGS="-fPIC" \
|
||||
-D CMAKE_CXX_FLAGS="-fPIC" \
|
||||
-D MPIEXEC=mpirun \
|
||||
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
|
||||
-D CMAKE_BUILD_TYPE:STRING=Release \
|
||||
-D USE_CUDA=1 \
|
||||
-D CMAKE_CUDA_FLAGS="-arch sm_80" \
|
||||
-D CMAKE_CUDA_HOST_COMPILER="$GCC_BIN/gcc" \
|
||||
-D USE_HDF5=1 \
|
||||
-D HDF5_DIRECTORY=$HDF5_DIR \
|
||||
-D USE_SILO=0 \
|
||||
-D SILO_DIRECTORY=$SILO_DIR \
|
||||
-D USE_TIMER=0 \
|
||||
~/LBPM-WIA
|
||||
|
||||
make VERBOSE=1 -j4 && make install
|
||||
|
||||
|
||||
|
|
@ -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 )
|
||||
|
|
|
|||
|
|
@ -107,8 +107,8 @@ int main (int argc, char *argv[])
|
|||
printf("Area ws = %f, Analytical = %f \n", Averages->aws, 4*PI*RADIUS*HEIGHT);
|
||||
printf("Area s = %f, Analytical = %f \n", Averages->As, 2*PI*RADIUS*(Nz-2));
|
||||
printf("Length wns = %f, Analytical = %f \n", Averages->lwns, 4*PI*RADIUS);
|
||||
printf("Geodesic curvature (wns) = %f, Analytical = %f \n", Averages->KGwns_global, 0.0);
|
||||
printf("Normal curvature (wns) = %f, Analytical = %f \n", Averages->KNwns_global, 1.0/RADIUS);
|
||||
printf("Geodesic curvature (wn) = %f, Analytical = %f \n", Averages->KGwns_global, 0.0);
|
||||
printf("Geodesic curvature (ws) = %f, Analytical = %f \n", Averages->KNwns_global, 4*PI);
|
||||
// printf("Cos(theta_wns) = %f, Analytical = %f \n",efawns/lwns,1.0*RADIUS/CAPRAD);
|
||||
printf("Interface Velocity = %f,%f,%f \n",Averages->vawn_global(0),Averages->vawn_global(1),Averages->vawn_global(2));
|
||||
printf("Common Curve Velocity = %f,%f,%f \n",Averages->vawns_global(0),Averages->vawns_global(1),Averages->vawns_global(2));
|
||||
|
|
|
|||
|
|
@ -74,7 +74,7 @@ int main(int argc, char **argv)
|
|||
while (timestep < Study.timestepMax && error > Study.tolerance){
|
||||
|
||||
timestep++;
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,false,0);//solve Poisson equtaion to get steady-state electrical potental
|
||||
IonModel.Run(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField); //solve for ion transport and electric potential
|
||||
|
||||
timestep++;//AA operations
|
||||
|
|
|
|||
|
|
@ -94,7 +94,7 @@ int main(int argc, char **argv)
|
|||
while (timestep < Study.timestepMax && error > Study.tolerance){
|
||||
|
||||
timestep++;
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,0);//solve Poisson equtaion to get steady-state electrical potental
|
||||
StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
|
||||
IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
|
||||
|
||||
|
|
|
|||
|
|
@ -69,7 +69,7 @@ int main(int argc, char **argv)
|
|||
int timeSave = int(PoissonSolver.TestPeriodicSaveInterval/PoissonSolver.TestPeriodicTimeConv);
|
||||
while (timestep<timeMax){
|
||||
timestep++;
|
||||
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,timestep);
|
||||
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,timestep);
|
||||
if (timestep%timeSave==0){
|
||||
if (rank==0) printf(" Time = %.3g[s]; saving electric potential and field\n",timestep*PoissonSolver.TestPeriodicTimeConv);
|
||||
PoissonSolver.getElectricPotential_debug(timestep);
|
||||
|
|
@ -78,7 +78,7 @@ int main(int argc, char **argv)
|
|||
}
|
||||
}
|
||||
else {
|
||||
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,1);
|
||||
PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,1);
|
||||
PoissonSolver.getElectricPotential_debug(1);
|
||||
PoissonSolver.getElectricField_debug(1);
|
||||
}
|
||||
|
|
|
|||
239
tests/lbpm_TwoPhase_analysis.cpp
Normal file
239
tests/lbpm_TwoPhase_analysis.cpp
Normal file
|
|
@ -0,0 +1,239 @@
|
|||
/*
|
||||
* 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);
|
||||
|
||||
/* move the solid surface by a voxel to improve contact angle measurement
|
||||
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;
|
||||
Averages->SDs(i,j,k) -= 1.0;
|
||||
}
|
||||
}
|
||||
}*/
|
||||
/* * * * * * * * * * * * * * * * * * * * * */
|
||||
// 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;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -32,10 +32,10 @@ int main( int argc, char **argv )
|
|||
int nprocs = comm.getSize();
|
||||
std::string SimulationMode = "production";
|
||||
// Load the input database
|
||||
auto db = std::make_shared<Database>( argv[1] );
|
||||
if (argc > 2) {
|
||||
SimulationMode = "legacy";
|
||||
}
|
||||
//auto db = std::make_shared<Database>( argv[1] );
|
||||
//if (argc > 2) {
|
||||
// SimulationMode = "legacy";
|
||||
//}
|
||||
|
||||
if ( rank == 0 ) {
|
||||
printf( "********************************************************\n" );
|
||||
|
|
@ -187,12 +187,13 @@ int main( int argc, char **argv )
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
PROFILE_STOP( "Main" );
|
||||
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
|
||||
auto level = db->getWithDefault<int>( "TimerLevel", 1 );
|
||||
NULL_USE(level);
|
||||
PROFILE_SAVE( file, level );
|
||||
*/
|
||||
// ****************************************************
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -94,7 +94,7 @@ int main(int argc, char **argv)
|
|||
while (timestep < Study.timestepMax){
|
||||
|
||||
timestep++;
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,timestep);//solve Poisson equtaion to get steady-state electrical potental
|
||||
PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
|
||||
StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
|
||||
IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
|
||||
|
||||
|
|
|
|||
|
|
@ -53,6 +53,7 @@ int main(int argc, char **argv)
|
|||
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
|
||||
SW = domain_db->getScalar<double>("Sw");
|
||||
auto READFILE = domain_db->getScalar<std::string>( "Filename" );
|
||||
auto MORPH_RADIUS = domain_db->getWithDefault<double>( "MorphRadius", 100000.0);
|
||||
|
||||
// Generate the NWP configuration
|
||||
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
|
||||
|
|
@ -122,7 +123,7 @@ int main(int argc, char **argv)
|
|||
comm.barrier();
|
||||
|
||||
// Run the morphological opening
|
||||
MorphDrain(SignDist, id, Dm, SW);
|
||||
MorphDrain(SignDist, id, Dm, SW, MORPH_RADIUS);
|
||||
|
||||
// calculate distance to non-wetting fluid
|
||||
if (domain_db->keyExists( "HistoryLabels" )){
|
||||
|
|
|
|||
|
|
@ -81,7 +81,7 @@ int main(int argc, char **argv)
|
|||
id = new signed char [N];
|
||||
Mask->Decomp(READFILE);
|
||||
Mask->CommInit();
|
||||
|
||||
|
||||
// Generate the NWP configuration
|
||||
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
|
||||
if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
|
||||
|
|
|
|||
|
|
@ -94,10 +94,11 @@ int main(int argc, char **argv)
|
|||
printf(" Measure the opening \n");
|
||||
sphere.MeasureObject();
|
||||
//sphere.ComputeScalar(Distance,0.f);
|
||||
printf(" Volume = %f (analytical = %f) \n", sphere.Vi,0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
|
||||
double error = fabs(sphere.Vi - 0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)))/ (0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
|
||||
/* Note 0.856 = (0.95)^3 */
|
||||
printf(" Volume = %f (analytical = %f) \n", sphere.Vi,0.256*0.33333333333333*0.856*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
|
||||
double error = fabs(sphere.Vi - 0.256*0.856*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)))/ (0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
|
||||
printf(" Relative error %f \n", error);
|
||||
if (error > 0.03){
|
||||
if (error > 0.05){
|
||||
toReturn = 10;
|
||||
printf("ERROR (test_morph): difference from analytical volume is too large\n");
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user