Repository: PMEAL/OpenPNM
Pull request #1148 (2018-11-14 23:52:37)

Unify meniscus models
Purcell and Sinusoidal are now wrappers to a general toroidal model which can also be used to model elliptical fibers. Purcell is valid for circular fibers which use elliptical model with scale_a = scale_b = r_toroid


Patch:
diff --git a/openpnm/models/physics/capillary_pressure.py b/openpnm/models/physics/capillary_pressure.py
index 2217987f4..630a2c25e 100644
--- a/openpnm/models/physics/capillary_pressure.py
+++ b/openpnm/models/physics/capillary_pressure.py
@@ -276,8 +276,9 @@ def ransohoff_snap_off(target,
     return value
 
 
-def purcell_bidirectional(target, r_toroid,
-                          num_points=1e2,
+def purcell_bidirectional(target,
+                          r_toroid=5e-6,
+                          num_points=1e3,
                           surface_tension='pore.surface_tension',
                           contact_angle='pore.contact_angle',
                           throat_diameter='throat.diameter',
@@ -326,7 +327,7 @@ def purcell_bidirectional(target, r_toroid,
         network['throat.temp_diameter'] = network[pore_diameter][conns[:, p]]
         key = 'throat.touch_pore_'+str(p)
         target.add_model(propname=key,
-                         model=pm.meniscus.toroidal,
+                         model=pm.meniscus.purcell,
                          mode='touch',
                          r_toroid=r_toroid,
                          num_points=num_points,
@@ -341,12 +342,11 @@ def purcell_bidirectional(target, r_toroid,
 
 
 def sinusoidal_bidirectional(target,
-                             num_points=1e2,
+                             r_toroid=5e-6,
+                             num_points=1e3,
                              surface_tension='pore.surface_tension',
                              contact_angle='pore.contact_angle',
                              throat_diameter='throat.diameter',
-                             throat_amplitude='throat.amplitude',
-                             throat_length='throat.length',
                              pore_diameter='pore.diameter'):
     r"""
     Computes the throat capillary entry pressure assuming the throat has a
@@ -364,6 +364,9 @@ def sinusoidal_bidirectional(target,
         controls the length of the calculated array, and also provides
         access to other necessary thermofluid properties
 
+    r_toroid : float or array_like
+        The radius of the toroid surrounding the pore
+
     num_points : float (Default 100)
         The number of divisions to make along the profile length to assess the
         meniscus properties in order to find the touch length.
@@ -379,13 +382,6 @@ def sinusoidal_bidirectional(target,
     throat_diameter : dict key (string)
         The dictionary key containing the throat diameter values to be used.
 
-    throat_amplitude : dict key (string)
-        The dictionary key containing the amplitude of variation in the throat
-        diameter about the mean.
-
-    throat_length : dict key (string)
-        The dictionary key containing the throat length values to be used.
-
     pore_diameter : dict key (string)
         The dictionary key containing the pore diameter values to be used.
     Notes
@@ -399,12 +395,11 @@ def sinusoidal_bidirectional(target,
         target.add_model(propname=key,
                          model=pm.meniscus.sinusoidal,
                          mode='touch',
+                         r_toroid=r_toroid,
                          num_points=num_points,
                          surface_tension=surface_tension,
                          contact_angle=contact_angle,
                          throat_diameter=throat_diameter,
-                         throat_amplitude=throat_amplitude,
-                         throat_length=throat_length,
                          touch_length='throat.temp_diameter')
         values[p] = target[key]
         target.remove_model(key)
diff --git a/openpnm/models/physics/meniscus.py b/openpnm/models/physics/meniscus.py
index 193e91472..c38a1b816 100644
--- a/openpnm/models/physics/meniscus.py
+++ b/openpnm/models/physics/meniscus.py
@@ -5,200 +5,107 @@ from openpnm.models.physics.capillary_pressure import _get_key_props
 logger = logging.getLogger(__name__)
 
 
-def toroidal(target,
-             mode='max',
-             r_toroid=5e-6,
-             target_Pc=None,
-             num_points=1e2,
-             surface_tension='pore.surface_tension',
-             contact_angle='pore.contact_angle',
-             throat_diameter='throat.diameter',
-             touch_length='throat.touch_length'):
-    r"""
-    Calculate the filling angle (alpha) for a given capillary pressure
-
-    Parameters
-    ----------
-    target : OpenPNM Object
-        The object for which these values are being calculated.  This
-        controls the length of the calculated array, and also provides
-        access to other necessary thermofluid properties.
-
-    mode : string (Default is 'max')
-        Determines what information to send back. Options are:
-        'max' : the maximum capillary pressure along the throat axis
-        'touch' : the maximum capillary pressure a meniscus can sustain before
-                  touching a solid feature
-        'men' : return the meniscus info for a target pressure
-
-    r_toroid : float or array_like
-        The radius of the toroid surrounding the pore
-
-    target_Pc : float
-        The target capillary pressure
-
-    num_points : float (Default 100)
-        The number of divisions to make along the profile length to assess the
-        meniscus properties in order to find target pressures, touch lengths,
-        minima and maxima.
-
-    surface_tension : dict key (string)
-        The dictionary key containing the surface tension values to be used. If
-        a pore property is given, it is interpolated to a throat list.
-
-    contact_angle : dict key (string)
-        The dictionary key containing the contact angle values to be used. If
-        a pore property is given, it is interpolated to a throat list.
-
-    throat_diameter : dict key (string)
-        The dictionary key containing the throat diameter values to be used.
-
-    touch_length : dict key (string)
-        The dictionary key containing the maximum length that a meniscus can
-        protrude into the connecting pore before touching a solid feature and
-        therfore invading
+def purcell(target,
+            mode='max',
+            target_Pc=None,
+            num_points=1e3,
+            r_toroid=5e-6,
+            throat_diameter='throat.diameter',
+            touch_length='throat.touch_length',
+            surface_tension='pore.surface_tension',
+            contact_angle='pore.contact_angle'):
+    r'''
+    Wrapper for the general toroidal model to implement the Purcell equation
+    for a torus with cylindrical profile and toroidal radius r_toroid.
 
     Notes
     -----
     This approach accounts for the converging-diverging nature of many throat
-    types.  Advancing the meniscus beyond the apex of the toroid requires an
+    types. Advancing the meniscus beyond the apex of the toroid requires an
     increase in capillary pressure beyond that for a cylindical tube of the
     same radius. The details of this equation are described by Mason and
     Morrow [1]_, and explored by Gostick [2]_ in the context of a pore network
     model.
-
-    References
-    ----------
-
-    .. [1] G. Mason, N. R. Morrow, Effect of contact angle on capillary
-           displacement curvatures in pore throats formed by spheres. J.
-           Colloid Interface Sci. 168, 130 (1994).
-    .. [2] J. Gostick, Random pore network modeling of fibrous PEMFC gas
-           diffusion media using Voronoi and Delaunay tessellations. J.
-           Electrochem. Soc. 160, F731 (2013).
-    """
-    network = target.project.network
-    phase = target.project.find_phase(target)
-    element, sigma, theta = _get_key_props(phase=phase,
-                                           diameter=throat_diameter,
-                                           surface_tension=surface_tension,
-                                           contact_angle=contact_angle)
-
-    x, R, rt, s, t = syp.symbols('x, R, rt, s, t')
-    # Equation of circle re-arranged for y
-    y = R*syp.sqrt(1 - (x/R)**2)
-    # Throat radius profile
-    r = rt + (R-y)
-    # Derivative of profile
-    rprime = r.diff(x)
-    # Filling angle
-    alpha = syp.atan(rprime)
-    # Radius of curvature of meniscus
-    rm = r/syp.cos(alpha+t)
-    # distance from center of curvature to meniscus contact point (Pythagoras)
-    a = syp.sqrt(rm**2 - r**2)
-    # angle between throat axis, meniscus center and meniscus contact point
-    gamma = syp.atan(r/a)
-    # Capillary Pressure
-    f = -2*s*syp.cos(alpha+t)/r
-    # Callable Functions
-    rx = syp.lambdify((x, R, rt), r, 'numpy')
-    fill_angle = syp.lambdify((x, R, rt), alpha, 'numpy')
-    Pc = syp.lambdify((x, R, rt, s, t), f, 'numpy')
-    rad_curve = syp.lambdify((x, R, rt, s, t), rm, 'numpy')
-    c2x = syp.lambdify((x, R, rt, s, t), a, 'numpy')
-    cap_angle = syp.lambdify((x, R, rt, s, t), gamma, 'numpy')
-
-    # Contact Angle
-    theta = np.deg2rad(theta)
-    # Network properties
-    throatRad = network[throat_diameter]/2
-    pos = np.arange(-r_toroid*0.999, r_toroid*0.999, r_toroid/num_points)
-    fiberRad = np.ones(len(throatRad))*r_toroid
-    # Now find the positions of the menisci along each throat axis
-    Y, X = np.meshgrid(throatRad, pos)
-    t_Pc = Pc(X, fiberRad, Y, sigma, theta)
-    # Values of minima and maxima
-    Pc_min = np.min(t_Pc, axis=0)
-    Pc_max = np.max(t_Pc, axis=0)
-    # Arguments of minima and maxima
-    a_min = np.argmin(t_Pc, axis=0)
-    a_max = np.argmax(t_Pc, axis=0)
-    if mode == 'max':
-        return Pc_max
-    elif mode == 'touch':
-        all_rad = rad_curve(X, fiberRad, Y, sigma, theta)
-        all_c2x = c2x(X, fiberRad, Y, sigma, theta)
-        all_cen = X + np.sign(all_rad)*all_c2x
-        dist = all_cen + np.abs(all_rad)
-        # Only count lengths where meniscus bulges into pore
-        dist[all_rad > 0] = 0.0
-        touch_len = network[touch_length]
-        mask = dist > touch_len
-        arg_touch = np.argmax(mask, axis=0)
-        # Make sure we only count ones that happen before max pressure
-        # And above min pressure (which will be erroneous)
-        arg_in_range = (arg_touch < a_max) * (arg_touch > a_min)
-        arg_touch[~arg_in_range] = a_max[~arg_in_range]
-        x_touch = pos[arg_touch]
-        # Return the pressure at which a touch happens
-        Pc_touch = Pc(x_touch, fiberRad, throatRad, sigma, theta)
-        return Pc_touch
-    elif target_Pc is None:
-        logger.exception(msg='Please supply a target capillary pressure' +
-                         ' when mode is "men"')
-    if np.abs(target_Pc) < 1.0:
-        target_Pc = 1.0
-
-    inds = np.indices(np.shape(t_Pc))
-    # Change values outside the range between minima and maxima to be those
-    # Values
-    mask = inds[0] < np.ones(len(pos))[:, np.newaxis]*a_min
-    t_Pc[mask] = (np.ones(len(pos))[:, np.newaxis]*Pc_min)[mask]
-    mask = inds[0] > np.ones(len(pos))[:, np.newaxis]*a_max
-    t_Pc[mask] = (np.ones(len(pos))[:, np.newaxis]*Pc_max)[mask]
-    # Find the argument at or above the target Pressure
-    mask = t_Pc >= target_Pc
-    arg_x = np.argmax(mask, axis=0)
-    # If outside range change to minima or maxima accordingly
-    arg_x[target_Pc < Pc_min] = a_min[target_Pc < Pc_min]
-    arg_x[target_Pc > Pc_max] = a_max[target_Pc > Pc_max]
-    xpos = pos[arg_x]
-    # Output
-    men_data = {}
-    men_data['pos'] = xpos
-    men_data['rx'] = rx(xpos, fiberRad, throatRad)
-    men_data['alpha'] = fill_angle(xpos, fiberRad, throatRad)
-    men_data['alpha_min'] = fill_angle(pos[a_min], fiberRad, throatRad)
-    men_data['alpha_max'] = fill_angle(pos[a_max], fiberRad, throatRad)
-    men_data['c2x'] = c2x(xpos, fiberRad, throatRad, sigma, theta)
-    men_data['gamma'] = cap_angle(xpos, fiberRad, throatRad, sigma, theta)
-    men_data['radius'] = rad_curve(xpos, fiberRad, throatRad, sigma, theta)
-    # xpos is relative to the throat center
-    men_data['center'] = (xpos + np.sign(men_data['radius'])*men_data['c2x'])
-    men_data['men_max'] = men_data['center'] - men_data['radius']
-    logger.info(mode+' calculated for Pc: '+str(target_Pc))
-    return men_data
+    '''
+    target['throat.scale_a'] = r_toroid
+    target['throat.scale_b'] = r_toroid
+    output = general_toroidal(target=target,
+                              mode=mode,
+                              profile_equation='elliptical',
+                              target_Pc=target_Pc,
+                              num_points=num_points,
+                              throat_diameter=throat_diameter,
+                              touch_length=touch_length,
+                              surface_tension=surface_tension,
+                              contact_angle=contact_angle)
+    return output
 
 
 def sinusoidal(target,
                mode='max',
                target_Pc=None,
-               num_points=1e2,
-               surface_tension='pore.surface_tension',
-               contact_angle='pore.contact_angle',
+               num_points=1e3,
+               r_toroid=5e-6,
                throat_diameter='throat.diameter',
-               throat_amplitude='throat.amplitude',
-               throat_length='throat.length',
+               pore_diameter='pore.diameter',
                touch_length='throat.touch_length',
-               **kwargs):
-    r"""
-    The profile of a throat is approximated with a sinusoidal function
-    that depends on the average of the connecting pore diameters and throat
-    diameter. It represents a converging-diverging geometry that has a minima
-    at the mid-point of the throat and produces similar behaviour to the
-    Purcell model but allows for a more slowly varying profile at the
-    ends of the throat.
+               surface_tension='pore.surface_tension',
+               contact_angle='pore.contact_angle'):
+    r'''
+    Wrapper for the general toroidal model to implement the a sinusoidal
+    profile. The quarter-wavelength is equal to toroidal radius r_toroid.
+    The amplitude is equal to the toroidal radius multiplied by the ratio of
+    the throat radius and average connecting pore radius.
+
+    Notes
+    -----
+    The capillary pressure equation for a sinusoidal throat is extended from
+    the Washburn equation as [1]_:
+
+    .. math::
+        P_c = -\frac{2\sigma(cos(\alpha + \theta))}{r(x)}
+
+    where alpha is:
+    .. math::
+        \alpha = arctan(\frac{dr}{dx})
+
+    References
+    ----------
+
+    .. [1] A. Forner-Cuenca et. al, Advanced Water Management in PEFCs:
+        Diffusion Layers with Patterned Wettability.
+        J. ECS. 163, 9, F1038-F1048 (2016).
+    '''
+    network = target.project.network
+    Dp_av = np.mean(network[pore_diameter][network['throat.conns']], axis=1)
+    scale_b = r_toroid*network[throat_diameter]/Dp_av
+    target['throat.scale_a'] = r_toroid
+    target['throat.scale_b'] = scale_b
+    output = general_toroidal(target=target,
+                              mode=mode,
+                              profile_equation='sinusoidal',
+                              target_Pc=target_Pc,
+                              num_points=num_points,
+                              throat_diameter=throat_diameter,
+                              touch_length=touch_length,
+                              surface_tension=surface_tension,
+                              contact_angle=contact_angle)
+    return output
+
+
+def general_toroidal(target,
+                     profile_equation='elliptical',
+                     mode='max',
+                     target_Pc=None,
+                     num_points=1e3,
+                     throat_scale_a='throat.scale_a',
+                     throat_scale_b='throat.scale_b',
+                     throat_diameter='throat.diameter',
+                     touch_length='throat.touch_length',
+                     surface_tension='pore.surface_tension',
+                     contact_angle='pore.contact_angle'):
+    r'''
+    The general model for meniscus properties inside a toroidal throat
 
     Parameters
     ----------
@@ -207,104 +114,105 @@ def sinusoidal(target,
         controls the length of the calculated array, and also provides
         access to other necessary thermofluid properties.
 
+    profile_equation : string (Default is 'elliptical')
+        options elliptical, sinusoidal.
+
     mode : string (Default is 'max')
         Determines what information to send back. Options are:
-        'max' : the maximum capillary pressure along the throat axis, does not
+        'max' : the maximum capillary pressure along the throat axis
         'touch' : the maximum capillary pressure a meniscus can sustain before
                   touching a solid feature
-        'men' : return the meniscus info for a target pressure
+        None : return the meniscus info for a target pressure
 
-    target_Pc : float (Default is None)
-        The target capillary pressure for use with mode 'men'
+    target_Pc : float
+        The target capillary pressure to return data for when mode is 'men'
 
     num_points : float (Default 100)
         The number of divisions to make along the profile length to assess the
         meniscus properties in order to find target pressures, touch lengths,
         minima and maxima.
 
-    surface_tension : dict key (string)
-        The dictionary key containing the surface tension values to be used. If
-        a pore property is given, it is interpolated to a throat list.
+    throat_scale_a : dict key (string)
+        The dictionary key containing the scale factor for adjusting the
+        profile along the throat axis (x).
 
-    contact_angle : dict key (string)
-        The dictionary key containing the contact angle values to be used. If
-        a pore property is given, it is interpolated to a throat list.
+    throat_scale_a : dict key (string)
+        The dictionary key containing the scale factor for adjusting the
+        profile perpendicular to the throat axis (y).
 
     throat_diameter : dict key (string)
-        The dictionary key containing the average throat diameter values.
-
-    throat_amplitude : dict key (string)
-        The dictionary key containing the amplitude of variation in the throat
-        diameter about the mean.
-
-    throat_length : dict key (string)
-        The dictionary key containing the throat length values to be used.
+        The dictionary key containing the throat diameter values to be used.
 
     touch_length : dict key (string)
         The dictionary key containing the maximum length that a meniscus can
         protrude into the connecting pore before touching a solid feature and
         therfore invading
 
-    Notes
-    -----
-    The capillary pressure equation for a sinusoidal throat is extended from
-    the Washburn equation as [1]_:
-
-    .. math::
-        P_c = -\frac{2\sigma(cos(\alpha + \theta))}{r(x)}
-
-    where alpha is:
-    .. math::
-        \alpha = arctan(\frac{dr}{dx})
-
-    References
-    ----------
+    surface_tension : dict key (string)
+        The dictionary key containing the surface tension values to be used. If
+        a pore property is given, it is interpolated to a throat list.
 
-    .. [1] A. Forner-Cuenca et. al, Advanced Water Management in PEFCs:
-        Diffusion Layers with Patterned Wettability.
-        J. ECS. 163, 9, F1038-F1048 (2016).
-    """
+    contact_angle : dict key (string)
+        The dictionary key containing the contact angle values to be used. If
+        a pore property is given, it is interpolated to a throat list.
+    '''
+    # Get data from dictionary keys
     network = target.project.network
     phase = target.project.find_phase(target)
-    element, sigma, theta = _get_key_props(phase=phase,
-                                           diameter=throat_diameter,
-                                           surface_tension=surface_tension,
-                                           contact_angle=contact_angle)
-    # Symbols
-    # position along throat, amplitude, throat rad, throat length, sigma, theta
-    x, A, rt, lt, s, t = syp.symbols('x, A, rt, lt, s, t')
-    # Equations
-    # Radius profile along throat length
-    y = A*syp.cos(2*syp.pi*x) + rt
-    # dr/dx used for filling angle
-    yprime = y.diff(x)/lt
+    (element,
+     surface_tension,
+     contact) = _get_key_props(phase=phase,
+                               diameter=throat_diameter,
+                               surface_tension=surface_tension,
+                               contact_angle=contact_angle)
+    # Contact Angle in radians
+    contact = np.deg2rad(contact)
+    # Network properties
+    throatRad = network[throat_diameter]/2
+    # Scaling parameters for throat profile
+    fa = target[throat_scale_a]
+    fb = target[throat_scale_b]
+    # Governing equations
+    x, a, b, rt, sigma, theta = syp.symbols('x, a, b, rt, sigma, theta')
+    if profile_equation == 'elliptical':
+        y = syp.sqrt(1 - (x/a)**2)*b
+    elif profile_equation == 'sinusoidal':
+        y = (syp.cos(syp.pi*x/(2*a)))*b
+    else:
+        logger.error('Profile equation is not valid')
+    # Throat radius profile
+    r = rt + (b-y)
+    # Derivative of profile
+    rprime = r.diff(x)
     # Filling angle
-    alpha = syp.atan(yprime)
-    # Meniscus Radius of curvature
-    R = y/syp.cos(alpha+t)
+    alpha = syp.atan(rprime)
+    # Radius of curvature of meniscus
+    rm = r/syp.cos(alpha+theta)
     # distance from center of curvature to meniscus contact point (Pythagoras)
-    a = syp.sqrt(R*R - y*y)
+    d = syp.sqrt(rm**2 - r**2)
     # angle between throat axis, meniscus center and meniscus contact point
-    gamma = syp.atan(y/a)
-    # Capillary Pressure function
-    f = -2*s*syp.cos(alpha+t)/y
-    # Callable expressions
-    rx = syp.lambdify((x, A, rt), y, 'numpy')
-    Pc = syp.lambdify((x, A, rt, lt, s, t), f, 'numpy')
-    rad_curve = syp.lambdify((x, A, rt, lt, s, t), R, 'numpy')
-    c2x = syp.lambdify((x, A, rt, lt, s, t), a, 'numpy')
-    fill_angle = syp.lambdify((x, A, rt, lt), alpha, 'numpy')
-    cap_angle = syp.lambdify((x, A, rt, lt, s, t), gamma, 'numpy')
-    # Contact Angle
-    theta = np.deg2rad(theta)
-    # Network properties
-    t_len = network[throat_length]
-    pos = np.arange(0.1, 0.9, 1/num_points)
-    r_amp = network[throat_amplitude]
-    r_ts = network[throat_diameter]/2
+    gamma = syp.atan(r/d)
+    # Capillary Pressure
+    p = -2*sigma*syp.cos(alpha+theta)/r
+    # Callable functions
+    rx = syp.lambdify((x, a, b, rt), r, 'numpy')
+    fill_angle = syp.lambdify((x, a, b, rt), alpha, 'numpy')
+    Pc = syp.lambdify((x, a, b, rt, sigma, theta), p, 'numpy')
+    rad_curve = syp.lambdify((x, a, b, rt, sigma, theta), rm, 'numpy')
+    c2x = syp.lambdify((x, a, b, rt, sigma, theta), d, 'numpy')
+    cap_angle = syp.lambdify((x, a, b, rt, sigma, theta), gamma, 'numpy')
+    # All relative positions along throat
+#    pos = np.arange(-0.999, 0.999, 1/num_points)
+    hp = int(num_points/2)
+    log_pos = np.logspace(-4, -1, hp+1)[:-1]
+    lin_pos = np.arange(0.1, 1.0, 1/hp)
+    half_pos = np.concatenate((log_pos, lin_pos))
+    pos = np.concatenate((-half_pos[::-1], half_pos))
     # Now find the positions of the menisci along each throat axis
-    Y, X = np.meshgrid(r_ts, pos)
-    t_Pc = Pc(X, r_amp, Y, t_len, sigma, theta)
+    Y, X = np.meshgrid(throatRad, pos)
+    X *= fa
+    # throat Capillary Pressure
+    t_Pc = Pc(X, fa, fb, Y, surface_tension, contact)
     # Values of minima and maxima
     Pc_min = np.min(t_Pc, axis=0)
     Pc_max = np.max(t_Pc, axis=0)
@@ -314,29 +222,30 @@ def sinusoidal(target,
     if mode == 'max':
         return Pc_max
     elif mode == 'touch':
-        all_rad = rad_curve(X, r_amp, Y, t_len, sigma, theta)
-        all_c2x = c2x(X, r_amp, Y, t_len, sigma, theta)
-        all_cen = X*t_len + np.sign(all_rad)*all_c2x
+        all_rad = rad_curve(X, fa, fb, Y, surface_tension, contact)
+        all_c2x = c2x(X, fa, fb, Y, surface_tension, contact)
+        all_cen = X + np.sign(all_rad)*all_c2x
         dist = all_cen + np.abs(all_rad)
         # Only count lengths where meniscus bulges into pore
         dist[all_rad > 0] = 0.0
-        touch_len = network[touch_length] + t_len/2
+        touch_len = network[touch_length]
         mask = dist > touch_len
         arg_touch = np.argmax(mask, axis=0)
         # Make sure we only count ones that happen before max pressure
         # And above min pressure (which will be erroneous)
         arg_in_range = (arg_touch < a_max) * (arg_touch > a_min)
         arg_touch[~arg_in_range] = a_max[~arg_in_range]
-        x_touch = pos[arg_touch]
+        x_touch = pos[arg_touch]*fa
         # Return the pressure at which a touch happens
-        Pc_touch = Pc(x_touch, r_amp, r_ts, t_len, sigma, theta)
+        Pc_touch = Pc(x_touch, fa, fb, throatRad, surface_tension, contact)
         return Pc_touch
     elif target_Pc is None:
         logger.exception(msg='Please supply a target capillary pressure' +
                          ' when mode is "men"')
     if np.abs(target_Pc) < 1.0:
         target_Pc = 1.0
-
+    # Find the position in-between the minima and maxima corresponding to
+    # the target pressure
     inds = np.indices(np.shape(t_Pc))
     # Change values outside the range between minima and maxima to be those
     # Values
@@ -350,19 +259,23 @@ def sinusoidal(target,
     # If outside range change to minima or maxima accordingly
     arg_x[target_Pc < Pc_min] = a_min[target_Pc < Pc_min]
     arg_x[target_Pc > Pc_max] = a_max[target_Pc > Pc_max]
-    xpos = pos[arg_x]
+    xpos = pos[arg_x]*fa
+    xmin = pos[a_min]*fa
+    xmax = pos[a_max]*fa
     # Output
     men_data = {}
     men_data['pos'] = xpos
-    men_data['rx'] = rx(xpos, r_amp, r_ts)
-    men_data['alpha'] = fill_angle(xpos, r_amp, r_ts, t_len)
-    men_data['c2x'] = c2x(xpos, r_amp, r_ts, t_len, sigma, theta)
-    men_data['gamma'] = cap_angle(xpos, r_amp, r_ts, t_len, sigma, theta)
-    men_data['radius'] = rad_curve(xpos, r_amp, r_ts, t_len, sigma, theta)
-    # xpos is relative to the start of the throat not the center
-    # Coop filling assumes center of throat
-    men_data['center'] = ((xpos-0.5)*t_len +
-                          np.sign(men_data['radius'])*men_data['c2x'])
+    men_data['rx'] = rx(xpos, fa, fb, throatRad)
+    men_data['alpha'] = fill_angle(xpos, fa, fb, throatRad)
+    men_data['alpha_min'] = fill_angle(xmin, fa, fb, throatRad)
+    men_data['alpha_max'] = fill_angle(xmax, fa, fb, throatRad)
+    men_data['c2x'] = c2x(xpos, fa, fb, throatRad, surface_tension, contact)
+    men_data['gamma'] = cap_angle(xpos, fa, fb, throatRad,
+                                  surface_tension, contact)
+    men_data['radius'] = rad_curve(xpos, fa, fb, throatRad,
+                                   surface_tension, contact)
+    # xpos is relative to the throat center
+    men_data['center'] = (xpos + np.sign(men_data['radius'])*men_data['c2x'])
     men_data['men_max'] = men_data['center'] - men_data['radius']
     logger.info(mode+' calculated for Pc: '+str(target_Pc))
     return men_data


Test patch:
diff --git a/tests/unit/algorithms/MixedPercolationTest.py b/tests/unit/algorithms/MixedPercolationTest.py
index e9c805c68..9632341db 100644
--- a/tests/unit/algorithms/MixedPercolationTest.py
+++ b/tests/unit/algorithms/MixedPercolationTest.py
@@ -437,14 +437,14 @@ class MixedPercolationTest:
                       model=op.models.geometry.throat_vector.pore_to_pore)
         water = op.phases.Water(network=pn)
         water['pore.contact_angle'] = 60
-        phys = op.physics.Standard(network=pn, phase=water, geometry=geo)
+        phys = op.physics.GenericPhysics(network=pn, phase=water, geometry=geo)
         r_tor = 5e-6
         phys.add_model(propname='throat.entry_pressure',
-                       model=op.models.physics.meniscus.toroidal,
+                       model=op.models.physics.meniscus.purcell,
                        r_toroid=r_tor,
                        mode='max')
         phys.add_model(propname='throat.meniscus',
-                       model=op.models.physics.meniscus.toroidal,
+                       model=op.models.physics.meniscus.purcell,
                        mode='men',
                        r_toroid=r_tor,
                        target_Pc=5000)
@@ -473,14 +473,14 @@ class MixedPercolationTest:
                       model=op.models.geometry.throat_vector.pore_to_pore)
         water = op.phases.Water(network=pn)
         water['pore.contact_angle'] = 100
-        phys = op.physics.Standard(network=pn, phase=water, geometry=geo)
+        phys = op.physics.GenericPhysics(network=pn, phase=water, geometry=geo)
         r_tor = 5e-6
         pmod = op.models.physics.capillary_pressure.purcell_bidirectional
         phys.add_model(propname='throat.entry_pressure',
                        model=pmod,
                        r_toroid=r_tor)
         phys.add_model(propname='throat.max_pressure',
-                       model=op.models.physics.meniscus.toroidal,
+                       model=op.models.physics.meniscus.purcell,
                        r_toroid=r_tor,
                        mode='max')
         phys['pore.entry_pressure'] = 0.0
diff --git a/tests/unit/models/physics/CapillaryPressureTest.py b/tests/unit/models/physics/CapillaryPressureTest.py
index 8b8ccf3f4..1db5c705e 100644
--- a/tests/unit/models/physics/CapillaryPressureTest.py
+++ b/tests/unit/models/physics/CapillaryPressureTest.py
@@ -91,9 +91,9 @@ class CapillaryPressureTest:
         f = op.models.physics.capillary_pressure.sinusoidal_bidirectional
         self.geo['pore.touch'] = (sp.random.random(self.geo.Np)+0.5)*0.1
         self.geo['throat.length'] = 1.0
-        self.geo['throat.amplitude'] = 0.25
         self.phys.add_model(propname='throat.bidirectional',
                             model=f,
+                            r_toroid=0.25,
                             surface_tension='pore.surface_tension',
                             contact_angle='pore.contact_angle',
                             pore_diameter='pore.touch')
diff --git a/tests/unit/models/physics/MeniscusTest.py b/tests/unit/models/physics/MeniscusTest.py
index e5d71fe32..f72c8477a 100644
--- a/tests/unit/models/physics/MeniscusTest.py
+++ b/tests/unit/models/physics/MeniscusTest.py
@@ -21,11 +21,11 @@ class MeniscusTest:
         r_tor = 1e-6
         self.geo['throat.touch_length'] = 2e-6
         phys.add_model(propname='throat.tor_max',
-                       model=pm.meniscus.toroidal,
+                       model=pm.meniscus.purcell,
                        mode='max',
                        r_toroid=r_tor)
         phys.add_model(propname='throat.tor_touch',
-                       model=pm.meniscus.toroidal,
+                       model=pm.meniscus.purcell,
                        mode='touch',
                        r_toroid=r_tor)
         assert sp.any(phys['throat.tor_touch'] < phys['throat.tor_max'])
@@ -69,12 +69,12 @@ class MeniscusTest:
                        model=pm.capillary_pressure.purcell,
                        r_toroid=r_tor)
         phys.add_model(propname='throat.tor_pressure',
-                       model=pm.meniscus.toroidal,
+                       model=pm.meniscus.purcell,
                        mode='max',
                        r_toroid=r_tor,
                        num_points=1000)
         phys.add_model(propname='throat.tor_meniscus',
-                       model=pm.meniscus.toroidal,
+                       model=pm.meniscus.purcell,
                        mode='men',
                        r_toroid=r_tor,
                        target_Pc=5000)
@@ -86,6 +86,26 @@ class MeniscusTest:
             if len(check) > 0:
                 assert 1 == 2
 
+    def test_general_toroidal(self):
+        phys = self.phys
+        r_tor = 1e-6
+        phys.add_model(propname='throat.purcell_pressure',
+                       model=pm.capillary_pressure.purcell,
+                       r_toroid=r_tor)
+        phys['throat.scale_a'] = r_tor
+        phys['throat.scale_b'] = r_tor
+        phys.add_model(propname='throat.general_pressure',
+                       model=pm.meniscus.general_toroidal,
+                       mode='max',
+                       num_points=1000)
+        a = sp.around(phys['throat.purcell_pressure'], 10)
+        b = sp.around(phys['throat.general_pressure'], 10)
+        assert sp.allclose(a, b)
+        h = phys.check_data_health()
+        for check in h.values():
+            if len(check) > 0:
+                assert 1 == 2
+
 
 if __name__ == '__main__':