From 9a6422c3630f51f37346d7e48f78dbae57bd39d2 Mon Sep 17 00:00:00 2001 From: Timothy Weigand <68024672+tmweigand@users.noreply.github.com> Date: Thu, 22 May 2025 22:29:27 -0400 Subject: [PATCH 1/3] updates drainage to optionally not operate on a multiphase.img --- src/pmmoto/domain_generation/multiphase.py | 13 +++++++++---- src/pmmoto/filters/equilibrium_distribution.py | 16 +++++++++++----- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/pmmoto/domain_generation/multiphase.py b/src/pmmoto/domain_generation/multiphase.py index c5db4d19..2cee505e 100644 --- a/src/pmmoto/domain_generation/multiphase.py +++ b/src/pmmoto/domain_generation/multiphase.py @@ -24,7 +24,7 @@ def update_img(self, img): """ self.img = img - def get_volume_fraction(self, phase: int) -> float: + def get_volume_fraction(self, phase: int, img=None) -> float: """ Calculate the volume fraction of a given phase in a multiphase image. @@ -34,7 +34,10 @@ def get_volume_fraction(self, phase: int) -> float: Returns: float: The volume fraction of the specified phase. """ - local_img = utils.own_img(self.subdomain, self.img) + if img is None: + img = self.img + + local_img = utils.own_img(self.subdomain, img) local_voxel_count = np.count_nonzero(local_img == phase) total_voxel_count = ( @@ -46,11 +49,13 @@ def get_volume_fraction(self, phase: int) -> float: total_voxels = np.prod(self.subdomain.domain.voxels) return total_voxel_count / total_voxels - def get_saturation(self, phase): + def get_saturation(self, phase: int, img=None) -> float: """ Calculate the saturation of a multiphase image """ - return self.get_volume_fraction(phase) / self.porous_media.porosity + if img is None: + img = self.img + return self.get_volume_fraction(phase, img) / self.porous_media.porosity @staticmethod def get_probe_radius(pc, gamma=1, contact_angle=0): diff --git a/src/pmmoto/filters/equilibrium_distribution.py b/src/pmmoto/filters/equilibrium_distribution.py index 7421fde8..79d7a7c3 100755 --- a/src/pmmoto/filters/equilibrium_distribution.py +++ b/src/pmmoto/filters/equilibrium_distribution.py @@ -42,12 +42,14 @@ def drainage( if contact_angle != 0: raise ValueError("The standard approach requires a zero contact angle!") approach = _standard_method + update_img = True elif method == "contact_angle": if contact_angle == 0: logging.warning( "The contact angle is zero. This will yield same results as the standard approach." ) approach = _contact_angle_method + update_img = False else: raise ValueError(f"{method} is not implemented. ") @@ -66,9 +68,13 @@ def drainage( ) # Update the phase distribution - multiphase.update_img( - np.where((morph == 1) & (w_connected == 2), 1, multiphase.img) - ) + if update_img: + multiphase.update_img( + np.where((morph == 1) & (w_connected == 2), 1, multiphase.img) + ) + mp_img = multiphase.img + else: + mp_img = np.where((morph == 1) & (w_connected == 2), 1, multiphase.img) if save: save_img_data_parallel( @@ -76,11 +82,11 @@ def drainage( ".", "_" ), multiphase.subdomain, - multiphase.img, + mp_img, ) # Store wetting phase saturation - w_saturation[n] = multiphase.get_saturation(2) + w_saturation[n] = multiphase.get_saturation(2, mp_img) if multiphase.subdomain.rank == 0: logging.info( From d2fc7b3de3dfab8d90dfcd4462c581cb6b6d7c9a Mon Sep 17 00:00:00 2001 From: Timothy Weigand <68024672+tmweigand@users.noreply.github.com> Date: Fri, 23 May 2025 10:58:06 -0400 Subject: [PATCH 2/3] add parameter to allow for elliptical inkbottle --- examples/drainage_inkbottle.py | 9 ++++- .../domain_generation/_domain_generation.pyx | 37 ++----------------- .../domain_generation/domain_generation.py | 7 ++-- 3 files changed, 15 insertions(+), 38 deletions(-) diff --git a/examples/drainage_inkbottle.py b/examples/drainage_inkbottle.py index a8194fe9..c389e684 100644 --- a/examples/drainage_inkbottle.py +++ b/examples/drainage_inkbottle.py @@ -82,7 +82,12 @@ def drain_ink_bottle(): reservoir_voxels=reservoir_voxels, ) - pm = pmmoto.domain_generation.gen_pm_inkbottle(sd) + # Scaling parameters for inkbottle + # Set to 1 for traditional inkbottle + r_y = 0.5 + r_z = 2 + + pm = pmmoto.domain_generation.gen_pm_inkbottle(sd, r_y, r_z) mp = pmmoto.domain_generation.gen_mp_constant(pm, 2) w_saturation_standard = pmmoto.filters.equilibrium_distribution.drainage( @@ -116,7 +121,7 @@ def drain_ink_bottle(): ) plt.xlabel("Wetting Phase Saturation") plt.ylabel("Capillary Pressure") - plt.savefig("examples/drainage_inkpottle.pdf") + plt.savefig("examples/drainage_inkbottle.pdf") plt.close() diff --git a/src/pmmoto/domain_generation/_domain_generation.pyx b/src/pmmoto/domain_generation/_domain_generation.pyx index d550dd52..e1e9df98 100755 --- a/src/pmmoto/domain_generation/_domain_generation.pyx +++ b/src/pmmoto/domain_generation/_domain_generation.pyx @@ -9,7 +9,6 @@ import numpy as np cimport numpy as cnp from numpy cimport uint8_t -from libcpp.unordered_map cimport unordered_map from libc.math cimport sin, cos cnp.import_array() @@ -24,7 +23,6 @@ __all__ = [ "gen_pm_sphere", "gen_pm_atom", "gen_inkbottle", - "gen_elliptical_inkbottle", ] def gen_pm_sphere(subdomain, spheres, kd: bool = False) -> np.ndarray: @@ -95,7 +93,7 @@ def gen_pm_atom(subdomain, atoms, kd: bool = False) -> np.ndarray: return gen_pm_sphere(subdomain, atoms) -def gen_inkbottle(double[:] x, double[:] y, double[:] z): +def gen_inkbottle(double[:] x, double[:] y, double[:] z, double r_y = 1.0, double r_z = 1.0): """ Generate pm for inkbottle test case. See Miller_Bruning_etal_2019 """ @@ -117,36 +115,9 @@ def gen_inkbottle(double[:] x, double[:] y, double[:] z): _grid[i,j,k] = 1 else: r = (0.01*cos(0.01*x[i]) + 0.5*sin(x[i]) + 0.75) - if (y[j]*y[j] + z[k]*z[k]) <= r*r: + ry = r*r_y + rz = r*r_z + if y[j]*y[j]/(ry*ry) + z[k]*z[k]/(rz*rz) <= 1: _grid[i,j,k] = 1 - return grid - - -def gen_elliptical_inkbottle(double[:] x, double[:] y, double[:] z): - """ - Generate ellipitical inkbottle test case. See Miller_Bruning_etal_2019 - """ - cdef int NX = x.shape[0] - cdef int NY = y.shape[0] - cdef int NZ = z.shape[0] - cdef int i, j, k - cdef double r - cdef double radiusY = 1.0 - cdef double radiusZ = 2.0 - - _grid = np.zeros((NX, NY, NZ), dtype=np.uint8) - cdef cnp.uint8_t [:,:,:] grid - - grid = _grid - - for i in range(0,NX): - for j in range(0,NY): - for k in range(0,NZ): - r = (0.01*cos(0.01*x[i]) + 0.5*sin(x[i]) + 0.75) - rY = r*radiusY - rz = r*radiusZ - if y[j]*y[j]/(rY*rY) + z[k]*z[k]/(rz*rz) <= 1: - _grid[i,j,k] = 1 - return grid \ No newline at end of file diff --git a/src/pmmoto/domain_generation/domain_generation.py b/src/pmmoto/domain_generation/domain_generation.py index 0e84cb43..107478aa 100644 --- a/src/pmmoto/domain_generation/domain_generation.py +++ b/src/pmmoto/domain_generation/domain_generation.py @@ -154,13 +154,14 @@ def gen_pm_atom_file( return pm -def gen_pm_inkbottle(subdomain): +def gen_pm_inkbottle(subdomain, r_y=1.0, r_z=1.0): """ - Generate an inkbottle with reservoirs + Generate an inkbottle with reservoirs. + To generate an elliptical inkbottle, r_y and r_z can be used. """ _img = _domain_generation.gen_inkbottle( - subdomain.coords[0], subdomain.coords[1], subdomain.coords[2] + subdomain.coords[0], subdomain.coords[1], subdomain.coords[2], r_y, r_z ) pm = porousmedia.gen_pm(subdomain, _img) utils.check_img_for_solid(subdomain, pm.img) From a1623c637571952148f25f6055293122e2040b61 Mon Sep 17 00:00:00 2001 From: Timothy Weigand <68024672+tmweigand@users.noreply.github.com> Date: Thu, 29 May 2025 21:27:55 -0400 Subject: [PATCH 3/3] ensure all drainage algorithm update the image progressively. fix bug in Schulz method and breaks into the two approaches --- examples/drainage_inkbottle.py | 23 +++++++- .../filters/equilibrium_distribution.py | 58 +++++++++++++++++-- 2 files changed, 75 insertions(+), 6 deletions(-) diff --git a/examples/drainage_inkbottle.py b/examples/drainage_inkbottle.py index c389e684..9f8c5602 100644 --- a/examples/drainage_inkbottle.py +++ b/examples/drainage_inkbottle.py @@ -65,7 +65,7 @@ def drain_ink_bottle(): voxels = (560, 120, 120) reservoir_voxels = 40 - subdomains = (8, 1, 1) + subdomains = (4, 1, 1) box = ((0.0, 14.0), (-1.5, 1.5), (-1.5, 1.5)) @@ -101,6 +101,20 @@ def drain_ink_bottle(): mp, capillary_pressure, gamma=1, contact_angle=20, method="contact_angle" ) + # Reset to fully wetted domain + mp = pmmoto.domain_generation.gen_mp_constant(pm, 2) + + # Bad method! + w_saturation_extended_contact_angle = ( + pmmoto.filters.equilibrium_distribution.drainage( + mp, + capillary_pressure, + gamma=1, + contact_angle=20, + method="extended_contact_angle", + ) + ) + # Save final state of multiphase image pmmoto.io.output.save_img_data_parallel( file_name="examples/drainage_inkbottle_img", @@ -119,8 +133,15 @@ def drain_ink_bottle(): ".", label="Contact Angle Method", ) + plt.plot( + w_saturation_extended_contact_angle, + capillary_pressure, + ".", + label="Extended Contact Angle Method", + ) plt.xlabel("Wetting Phase Saturation") plt.ylabel("Capillary Pressure") + plt.legend() plt.savefig("examples/drainage_inkbottle.pdf") plt.close() diff --git a/src/pmmoto/filters/equilibrium_distribution.py b/src/pmmoto/filters/equilibrium_distribution.py index 79d7a7c3..56bf80b0 100755 --- a/src/pmmoto/filters/equilibrium_distribution.py +++ b/src/pmmoto/filters/equilibrium_distribution.py @@ -16,7 +16,7 @@ def drainage( capillary_pressures, gamma=1, contact_angle=0, - method: Literal["standard", "contact_angle"] = "standard", + method: Literal["standard", "contact_angle", "extended_contact_angle"] = "standard", save=False, ): """ @@ -49,7 +49,14 @@ def drainage( "The contact angle is zero. This will yield same results as the standard approach." ) approach = _contact_angle_method - update_img = False + update_img = True + elif method == "extended_contact_angle": + if contact_angle == 0: + logging.warning( + "The contact angle is zero. This will yield same results as the standard approach." + ) + approach = _extended_contact_angle_method + update_img = True else: raise ValueError(f"{method} is not implemented. ") @@ -83,6 +90,7 @@ def drainage( ), multiphase.subdomain, mp_img, + additional_img={"morph": morph}, ) # Store wetting phase saturation @@ -101,6 +109,9 @@ def drainage( def _standard_method(multiphase, capillary_pressure, gamma, contact_angle): """ This method for drainage follows Hilpert and Miller 2001 + The radius(r) is defined as: + r = 2*gamma / p_c + where gamma is the surface tensions and p_c is the capillary pressure. """ # Compute morphological changes based on capillary pressure @@ -123,17 +134,54 @@ def _standard_method(multiphase, capillary_pressure, gamma, contact_angle): def _contact_angle_method(multiphase, capillary_pressure, gamma, contact_angle): + """ + This method for drainage follows Schulz and Becker 2007. + The radius(r) is defined as: + r = 2*gamma*cos(theta) / p_c + where gamma is the surface tensions, theta is the contact angle, + and p_c is the capillary pressure. + """ + + # Compute morphological changes based on capillary pressure + radius = multiphase.get_probe_radius(capillary_pressure, gamma, contact_angle) + + # Check if radius is larger than resolution + if radius < min(multiphase.subdomain.domain.resolution): + morph = np.zeros_like(multiphase.pm_img) + else: + morph = porosimetry( + subdomain=multiphase.subdomain, + porous_media=multiphase.porous_media, + radius=radius, + inlet=True, + multiphase=multiphase, + mode="hybrid", + ) + + return morph + + +def _extended_contact_angle_method( + multiphase, capillary_pressure, gamma, contact_angle +): """ This method for drainage follows Schulz and Wargo 2015. - This paper seems to have many typos/ errors? - ..For example Eqations 1 and 5 show dimeter but should be radii. + In this work the radius (r) of erosion step includes the contact angle + but the dilation does not. + + + This method does not yield good results """ # Compute morphological changes based on capillary pressure erosion_radius = multiphase.get_probe_radius( capillary_pressure, gamma, contact_angle ) - dilation_radius = multiphase.get_probe_radius(capillary_pressure, gamma) + + # Mutiply again by cosine of contact angle + dilation_radius = multiphase.get_probe_radius( + capillary_pressure, gamma, contact_angle + ) * np.abs(np.cos(np.deg2rad(contact_angle))) radius = [erosion_radius, dilation_radius]