From 66be77aeaed8b43db8d2fa3fc887deb394c58bd7 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 14 Dec 2021 16:46:00 -0500 Subject: [PATCH] update deficit curvature test --- analysis/TwoPhase.cpp | 5 +- docs/source/examples/analysis/TwoPhase.rst | 135 ++++++++++++++++++ docs/source/examples/color/fractionalFlow.rst | 125 ++++++++++++++++ docs/source/examples/running.rst | 3 +- tests/TestInterfaceSpeed.cpp | 4 +- 5 files changed, 266 insertions(+), 6 deletions(-) create mode 100644 docs/source/examples/analysis/TwoPhase.rst create mode 100644 docs/source/examples/color/fractionalFlow.rst diff --git a/analysis/TwoPhase.cpp b/analysis/TwoPhase.cpp index bfa644b8..39ccbe4d 100644 --- a/analysis/TwoPhase.cpp +++ b/analysis/TwoPhase.cpp @@ -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; diff --git a/docs/source/examples/analysis/TwoPhase.rst b/docs/source/examples/analysis/TwoPhase.rst new file mode 100644 index 00000000..31bead84 --- /dev/null +++ b/docs/source/examples/analysis/TwoPhase.rst @@ -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 + + + diff --git a/docs/source/examples/color/fractionalFlow.rst b/docs/source/examples/color/fractionalFlow.rst new file mode 100644 index 00000000..c22e521d --- /dev/null +++ b/docs/source/examples/color/fractionalFlow.rst @@ -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) + diff --git a/docs/source/examples/running.rst b/docs/source/examples/running.rst index 70a36dfa..30cc5801 100644 --- a/docs/source/examples/running.rst +++ b/docs/source/examples/running.rst @@ -18,4 +18,5 @@ a basic introduction to working with LBPM. morphology/* color/* - + + analysis/* diff --git a/tests/TestInterfaceSpeed.cpp b/tests/TestInterfaceSpeed.cpp index 4036a205..1e8e271c 100644 --- a/tests/TestInterfaceSpeed.cpp +++ b/tests/TestInterfaceSpeed.cpp @@ -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));