Repository: PMEAL/OpenPNM
Pull request #978 (2018-07-11 06:53:31)

Minor tweak to use 64 bit integers in A matrix
This was added to force scipy to use the 64 bit version of umfpack and hence increase the allowed memory.  It works! 


Patch:
diff --git a/openpnm/algorithms/GenericTransport.py b/openpnm/algorithms/GenericTransport.py
index b2f6da951..04b8d59d1 100644
--- a/openpnm/algorithms/GenericTransport.py
+++ b/openpnm/algorithms/GenericTransport.py
@@ -6,6 +6,7 @@ from scipy.spatial import cKDTree
 from openpnm.topotools import iscoplanar
 from openpnm.algorithms import GenericAlgorithm
 from openpnm.utils import logging
+import inspect
 # Check if petsc4py is available
 import importlib
 if (importlib.util.find_spec('petsc4py') is not None):
@@ -434,8 +435,19 @@ class GenericTransport(GenericAlgorithm):
             x = sls.solve(ls)
             del(ls)  # Clean
         else:
+            A = A.tocsr()
+            A.indices = A.indices.astype(np.int64)
+            A.indptr = A.indptr.astype(np.int64)
             solver = getattr(sprs.linalg, self.settings['solver'])
-            x = solver(A=A.tocsr(), b=b)
+            if 'tol' in inspect.getargspec(solver)[0]:
+                norm_A = sprs.linalg.norm(self._A)
+                norm_b = np.linalg.norm(self._b)
+                tol = min(norm_A, norm_b)*1e-06
+                x = solver(A=A, b=b, tol=tol)
+            else:
+                x = solver(A=A, b=b)
+        if type(x) == tuple:
+            x = x[0]
         return x
 
     def results(self):
diff --git a/openpnm/algorithms/ReactiveTransport.py b/openpnm/algorithms/ReactiveTransport.py
index 18eb7e1a1..16df5ca55 100644
--- a/openpnm/algorithms/ReactiveTransport.py
+++ b/openpnm/algorithms/ReactiveTransport.py
@@ -94,7 +94,7 @@ class ReactiveTransport(GenericTransport):
         self._update_physics()
 
         x = self._run_reactive(x=x)
-        return x
+        self[self.settings['quantity']] = x
 
     def _run_reactive(self, x):
         if x is None:
diff --git a/requirements.txt b/requirements.txt
index 5f6f2ae77..4db0361ac 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,3 +10,4 @@ h5py
 porespy
 transforms3d
 flatdict
+scikit-umfpack
diff --git a/setup.py b/setup.py
index 279a363fc..0b36b2695 100644
--- a/setup.py
+++ b/setup.py
@@ -54,7 +54,8 @@ setup(
         'h5py',
         'porespy',
         'transforms3d',
-        'flatdict'],
+        'flatdict',
+        'scikit-umfpack'],
 
     author='OpenPNM Team',
     author_email='jgostick@uwaterloo.ca',


Test patch:
diff --git a/tests/unit/algorithms/GenericTransportTest.py b/tests/unit/algorithms/GenericTransportTest.py
index 9aa8d5359..f211841d1 100644
--- a/tests/unit/algorithms/GenericTransportTest.py
+++ b/tests/unit/algorithms/GenericTransportTest.py
@@ -55,6 +55,19 @@ class GenericTransportTest:
         y = sp.unique(sp.around(alg['pore.mole_fraction'], decimals=3))
         assert sp.all(x == y)
 
+    def test_two_value_conditions_cg(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        alg.settings['conductance'] = 'throat.diffusive_conductance'
+        alg.settings['quantity'] = 'pore.mole_fraction'
+        alg.set_value_BC(pores=self.net.pores('top'), values=1)
+        alg.set_value_BC(pores=self.net.pores('bottom'), values=0)
+        alg.settings['solver'] = 'cg'
+        alg.run()
+        x = [0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
+        y = sp.unique(sp.around(alg['pore.mole_fraction'], decimals=3))
+        assert sp.all(x == y)
+
     def test_one_value_one_rate(self):
         alg = op.algorithms.GenericTransport(network=self.net,
                                              phase=self.phase)