import numpy as np
import openpnm as op
import matplotlib.pyplot as plt
%config InlineBackend.figure_formats = ['svg']
np.random.seed(10)
ws = op.Workspace()
ws.settings["loglevel"] = 40
np.set_printoptions(precision=5)

"""
Explanation: Mercury Intrusion Porosimetry
Simulating capillary pressure curves using Porosimetry
Start by importing OpenPNM.
End of explanation
"""


pn = op.network.Cubic(shape=[20, 20, 20], spacing=0.00005)

"""
Explanation: Next, create a simple cubic network with 20 pores per side and a spacing of 50 um
End of explanation
"""


geo = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)

"""
Explanation: The network object only contains topological and spatial information.  We need to assign some pore and throat sizes, which can be conveniently done by creating a SpheresAndCylinders geometry object.
End of explanation
"""


hg = op.phases.Mercury(network=pn, name='mercury')

"""
Explanation: The SpheresAndCylinders class assigns pores as random values between 0.2 and 0.7 of the lattice spacing, then finds all other geometric information assuming the pores are spheres and the throats are cylinders.
Next we need to create a phase object which contains the thermo-physical properties of the invading fluid, such as surface tension:
End of explanation
"""


phys = op.physics.GenericPhysics(network=pn, phase=hg, geometry=geo)
phys.add_model(propname='throat.entry_pressure',
               model=op.models.physics.capillary_pressure.washburn)

"""
Explanation: Lastly, we need to compute the capillary entry pressure of the throats in the network.  The OpenPNM models library has a few common capillary pressure models, including the Washburn equation:
$$ P_C = \frac{-2\sigma cos(\theta)}{R_T} $$
To use this model we should create a physics object, and use the add_model method as follows:
End of explanation
"""


print(pn.project.grid)

"""
Explanation: Note that we can inspect our project to see the interrelationships between all the object with:
End of explanation
"""


mip = op.algorithms.Porosimetry(network=pn, phase=hg)
mip.set_inlets(pores=pn.pores('front'))
mip.run()

"""
Explanation: The grid shows us that phys_01 is associated with the mercury phase, and geo_01.  This means that when calculating the throat entry pressure using the Washburn equation above, it will fetch the throat radius ($R_T$) from geo_01 and the surface tension and contact angle from mercury.
Now that all the required objects are setup, we can perform the capillary pressure curve simulation.  OpenPNM contains both InvasionPercolation and OrdinaryPercolation classes.  The porosimetry experiment is ordinary percolation with access limitations.  This means that a fixed pressure is applied to the invading fluid and all pores and throat that can be penetrated at that pressure are set as possibly invaded, then pores and throats are set to invaded only if they are physically connected to the source of invading fluid directly or though a path of already invading pores and throats.  Since this process is simulated very frequently, OpenPNM includes a class called Porosimetry that is a subclass of OrdinaryPercolation, with several useful methods added.  It's use is illustrated below:
End of explanation
"""


mip.plot_intrusion_curve()

Pc, Snwp = mip.get_intrusion_data()
print(Pc, Snwp)

"""
Explanation: The meaning of this warning message will be analyzed below, but first let's take a quick look at the result using the built-in plotting method:
End of explanation
"""


fig, ax = plt.subplots()
ax.semilogx(Pc, Snwp, 'r*-')
ax.set_xlabel("Capillary pressure (Pa)")
ax.set_ylabel("Saturation")

"""
Explanation: With the above data in the form of arrays it's possible to cut&paste into Excel, or to use a Python plotting package such as Matplotlib to make plots with your desired style:
End of explanation
"""


