import scipy as sp
import openpnm as op

"""
Explanation: Part 1: Import Networks from Statoil Files
This example explains how to use the OpenPNM.Utilies.IO.Statoil class to import a network produced by the Maximal Ball network extraction code developed by Martin Blunt's group at Imperial College London. The code is available from him upon request, but they offer a small library of pre-extracted networks on their website.
End of explanation
"""


from pathlib import Path
path = Path('../fixtures/ICL-Sandstone(Berea)/')
project = op.io.Statoil.load(path=path, prefix='Berea')
pn = project.network
pn.name = 'berea'

"""
Explanation: The following assumes that the folder containing the 'dat' files is in a directory called 'fixtures' in the same directory as this script.  You can also enter a full path to the files.
End of explanation
"""


print(pn)

"""
Explanation: This import class extracts all the information contained in the 'Statoil' files, such as sizes, locations and connectivity. Note that the io classes return a project object, and the network itself can be accessed using the network attribute.  The following printout display which information was contained in the file:
End of explanation
"""


project.export_data(network=pn, filename='imported_statoil')

"""
Explanation: At this point, the network can be visualized in Paraview. A suitable '.vtp' file can be created with:
End of explanation
"""


print('Number of pores before trimming: ', pn.Np)
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
print('Number of pores after trimming: ', pn.Np)

"""
Explanation: The resulting network is shown below:

Clean up network topology
Although it's not clear in the network image, there are a number of isolated and disconnected pores in the network.  These are either naturally part of the sandstone, or artifacts of the Maximal Ball algorithm.  In any event, these must be removed before proceeding since they cause problems for the matrix solver.  The easiest way to find these is to use the check_network_health function on the network object.  This will return a dictionary with several key attributes, including a list of which pores are isolated. These can then be trimmed using the trim function in the topotools module.
End of explanation
"""


water = op.phases.Water(network=pn)

"""
Explanation: Dealing with Inlet and Outlet Pores
When importing Statoil networks, OpenPNM must perform some 'optimizations' to make the network compatible.  The main problem is that the original network contains a large number of throats connecting actual internal pores to fictitious 'reservoir' pores.  OpenPNM strips away all these throats since 'headless throats' break the graph theory representation.  OpenPNM then labels the real internal pores as either 'inlet' or 'outlet' if they were connected to one of these fictitious reservoirs.  
It is fairly simple to add a new pores to each end of the domain and stitch tehm to the internal pores labelled 'inlet' and 'outlet', but this introduces a subsequent complication that the new pores don't have any geometry properties.  For this example, we will not add boundary pores, but just the pores on the inlet and outlet faces.  
Part 2: Calculating Permeability of the Network
Setup Geometry, Phase, and Physics Objects
In OpenPNM 2+ it is optional to define Geometry and Physics objects (These are really only necessary for simulations with diverse geometrical properties in different regions, resulting in different physical processes in each region, such as multiscale networks for instance).  It is still necessary to define Phase objects:
End of explanation
"""


pn['pore.diameter'] = 2*pn['pore.radius']
pn['throat.diameter'] = 2*pn['throat.radius']
pn.add_model(propname='pore.area',
             model=op.models.geometry.pore_area.sphere)
pn.add_model(propname='throat.area',
             model=op.models.geometry.throat_area.cylinder)
pn.add_model(propname='throat.endpoints',
             model=op.models.geometry.throat_endpoints.spherical_pores)
pn.add_model(propname='throat.conduit_lengths',
             model=op.models.geometry.throat_length.conduit_lengths)

"""
Explanation: Also, we must add a few geometrical properties that were not included in the extracted network:
End of explanation
"""


water.add_model(propname='throat.hydraulic_conductance',
                model=op.models.physics.hydraulic_conductance.hagen_poiseuille)

"""
Explanation: Apply Pore-Scale Models
We must add the hagen-poiseuille model for calculating the conductance.  In OpenPNM 2+ it is possible to add Physics models to Phase objects, which is often simpler than than applying the same model to multiple Physics.
End of explanation
"""


print(water['throat.hydraulic_conductance'])

"""
Explanation: Recall that boundary pores and throats had no geometrical properties associated with them, so the   hydraulic conductances of boundary throats will be undefined (filled with NaNs):
End of explanation
"""


flow = op.algorithms.StokesFlow(network=pn, phase=water)
flow.set_value_BC(pores=pn.pores('inlets'), values=200000)
flow.set_value_BC(pores=pn.pores('outlets'), values=100000)
flow.run()

"""
Explanation: Run StokesFlow Algorithm
Finally, we can create a StokesFlow object to run some fluid flow simulations:
End of explanation
"""


# Get the average value of the fluid viscosity
mu = sp.mean(water['pore.viscosity'])
# Specify a pressure difference (in Pa)
delta_P = 100000
# Using the rate method of the StokesFlow algorithm
Q = sp.absolute(flow.rate(pores=pn.pores('inlets')))
# Because we know the inlets and outlets are at x=0 and x=X
Lx = sp.amax(pn['pore.coords'][:, 0]) - sp.amin(pn['pore.coords'][:, 0])
A = Lx*Lx  # Since the network is cubic Lx = Ly = Lz
K = Q*mu*Lx/(delta_P*A)
print(K)

"""
Explanation: The resulting pressure field can be visualized in Paraview, giving the following:

Determination of Permeability Coefficient
The way to calculate K is the determine each of the values in Darcy's law manually and solve for K, such that $$ K = \frac{Q\mu L} {\Delta P A} $$
End of explanation
"""


