Update Schulz Method

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))
 
@@ -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(
@@ -96,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",
@@ -114,9 +133,16 @@ 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.savefig("examples/drainage_inkpottle.pdf")
+        plt.legend()
+        plt.savefig("examples/drainage_inkbottle.pdf")
         plt.close()
 
 

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

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)

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

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,
 ):
     """
@@ -42,12 +42,21 @@ 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 = 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. ")
 
@@ -66,21 +75,26 @@ 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(
                 f"drainage_results/capillary_pressure_{capillary_pressure:.3f}".replace(
                     ".", "_"
                 ),
                 multiphase.subdomain,
-                multiphase.img,
+                mp_img,
+                additional_img={"morph": morph},
             )
 
         # 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(
@@ -95,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
@@ -117,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]