update deficit curvature test

This commit is contained in:
James McClure 2021-12-14 16:46:00 -05:00
parent 012c1814c6
commit 66be77aeae
5 changed files with 266 additions and 6 deletions

View File

@ -1472,17 +1472,16 @@ void TwoPhase::Reduce() {
van_global(2) = van_global(2) / nwp_volume_global;
}
// Normalize surface averages by the interfacial area
/*
if (awn_global > 0.0) {
Jwn_global /= awn_global;
Kwn_global /= awn_global;
//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;

View File

@ -0,0 +1,135 @@
*************************
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.
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

View 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)

View File

@ -18,4 +18,5 @@ a basic introduction to working with LBPM.
morphology/*
color/*
analysis/*

View File

@ -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));