periodic flow simulation

gen/generate.py

@@ -7,8 +7,8 @@ from netflow import *
 from time import perf_counter
         
 def main():
-    if len(sys.argv) != 5:
-        raise AssertionError("Usage: python3 generate.py <nthreads> <n1> <n2> <n3>")
+    if len(sys.argv) != 6:
+        raise AssertionError("Usage: python3 generate.py <nthreads> <n1> <n2> <n3> <save network?: y/n>")
     nthreads = int(sys.argv[1])
     
     print("Using %d thread(s)" % nthreads)
@@ -26,7 +26,8 @@ def main():
     
     start = perf_counter()
     dendro = netflow.generate_dendrogram(basenet=basenet, targetsize=target, \
-                                        cutoff=cutoff, sd=42, nthreads=nthreads, mute=False)
+                                        cutoff=cutoff, sd=42, nthreads=nthreads, \
+                                        maxIters=6, throatTol=0.01, mute=False)
     end = perf_counter()
     elapsed = end - start
     print("Elapsed time: %e s" % elapsed)
@@ -35,9 +36,11 @@ def main():
     _ = netgen.plot_network(dendro)
     plt.savefig('plots/dendro.png')
     """
-    """
-    netflow.save_network_to('../visualize/network/dendro.h5', dendro)
-    """
+    if (sys.argv[5].lower() == 'y'):
+        path = 'networks/dendro_{}_{}_{}.h5'.format(*target)
+        netflow.save_network_to(path, dendro)
+        print("network saved to: {}".format(path))
+        
     print("Generated network statistics:")
     print("Throats: %d" % len(dendro.throats))
     print("Pores:   %d" % len(dendro.pores))

netflow/netgen.py

@@ -635,7 +635,7 @@ def cell_list_scaling(basenet: Network, targetsize: List[int], \
     
 def generate_dendrogram(basenet: Network, targetsize: List[int], \
     cutoff: float = float('inf'), sd: int = None, nthreads: int = mp.cpu_count(), \
-    mute: bool = False) -> Network:
+    maxIters: int = 6, throatTol: float = 0.01, mute: bool = False) -> Network:
     """Generate a new spatially periodic network.
 
     Based on an existing network basenet, generate a new network of size
@@ -776,19 +776,15 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
     
     if (not mute): print("Total throats: %d" % totalThroats)
     
-    # Additional parameters:
-    maxIters = 6  # Maximum number of iterations until we give up
-    throatTolerance = 0.01  # Max. 1% may not be realized
-    
     iterCount = 0
     throatsLeft = totalThroats
-    throatThreshold = int(throatTolerance * totalThroats)
+    throatThreshold = int(throatTol * totalThroats)
     
     # Pores in interior that are not fully-connected yet
     poresRemain = pores[:n]
     while (iterCount < maxIters and throatsLeft > throatThreshold):
         
-        if (not mute): 
+        if (not mute):
             print("Throats left: %d (%.1f%%)" % \
                  (throatsLeft, throatsLeft / totalThroats * 100.))
         workers = []

petsc-solve/cwrapper/wrapper.pyx

@@ -30,7 +30,7 @@ def __finalize_petsc():
 # A: scipy.sparse.csr_matrix
 # b: numpy.ndarray
 def solve_py(A, b, PetscScalar rtol, PetscInt maxiter, PetscInt n_rows):
-    __init_petsc()
+    ierr = __init_petsc()
     
     # Output: solution vector
     cdef PetscScalar [::1] solution = np.empty(n_rows, dtype=np.dtype("d"))
@@ -43,11 +43,11 @@ def solve_py(A, b, PetscScalar rtol, PetscInt maxiter, PetscInt n_rows):
     cdef const PetscInt [::1] row_ptr = A.indptr
     cdef const PetscInt [::1] col_indices = A.indices
     
-    # TODO: Check error code returned by solve 
+    # TODO: Check error code returned by solve
     solve(&data[0], &col_indices[0], &row_ptr[0], &rhs[0], &solution[0], \
                   rtol, maxiter, n_rows)
 	
-    __finalize_petsc()
+    ierr = __finalize_petsc()
     # Convert memory view into np.ndarray and return
     return np.asarray(solution)
 	

petsc-solve/solve_flow.py

+from mpi4py import MPI
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+root = 0
+
+import sys
+sys.path.append('../')
+
+from netflow import *
+
+def main():
+    
+    if (len(sys.argv) != 2):
+        if (rank == root):
+            print("Usage: python3 solve_flow.py <path-to-network>")
+        comm.Abort(1)
+    
+    net = None
+    
+    # Only initialize on root rank
+    if (rank == root):
+        netPath = sys.argv[1]
+        try:
+            net = netflow.load_network_from(netPath, isBaseNet=False)
+        except Exception as e:
+            print("Error: {}".format(e))
+            comm.Abort(1)
+        
+        print("Pores: %d" % (len(net.pores)))
+        
+    start = MPI.Wtime()
+    pA, QA = netflow.solve_flow_periodic(network=net, \
+                solver=netflow.Solver.PETSC, mu=1e-3, c=1, P2L=-1e3)
+    end = MPI.Wtime()
+    
+    if rank == root:
+        print("PETSc time: %e s" % (end - start))
+    
+    start = MPI.Wtime()
+    pB, QB = netflow.solve_flow_periodic(network=net, \
+                solver=netflow.Solver.CG, mu=1e-3, c=1, P2L=-1e3)
+    end = MPI.Wtime()
+    
+    if rank == root:
+        print("Other time: %e s" % (end - start))
+        
+        diff = 0.
+        count = 0
+        for pore in net.pores:
+            if pore in pA:
+                count += 1
+                diff += (pA[pore] - pB[pore])**2
+        print("MSE: %e" % (diff / count))
+        
+if __name__ == '__main__':
+    main()

petsc-solve/src/solve.c

@@ -24,7 +24,7 @@ PetscErrorCode solve(const PetscScalar *data, const PetscInt *col_indices,
                PetscInt maxiter, PetscInt n_rows) {
 	
     PetscErrorCode ierr;
-    
+        
     const PetscMPIInt root = 0;
     // Fetch number of processes and local rank
     PetscMPIInt size, rank;

petsc-solve/test.py

@@ -42,7 +42,7 @@ def main():
     
     # Initialization on root only
     if (rank == root):
-        basenet = netflow.load_network_from('../netflow/network/network.h5')
+        basenet = netflow.load_network_from('../netflow/network/network.h5', isBaseNet=False)
         for pore in basenet.pores:
             if pore.label[:2] == 'in':
                 inpores.add(pore)
@@ -64,7 +64,7 @@ def main():
         if rank == root:
             print("Solver: %s Time: %e s" % (solver.name, end - start))
             
-            plot_flux(inpores, outpores, p, Q, solver.name)
+            #plot_flux(inpores, outpores, p, Q, solver.name)
         comm.Barrier()
 	
 if __name__ == '__main__':

thesis/thesis.log

-This is pdfTeX, Version 3.14159265-2.6-1.40.20 (TeX Live 2019/Debian) (preloaded format=pdflatex 2021.4.27)  30 NOV 2021 22:10
+This is pdfTeX, Version 3.14159265-2.6-1.40.20 (TeX Live 2019/Debian) (preloaded format=pdflatex 2021.4.27)  1 DEC 2021 12:30
 entering extended mode
  restricted \write18 enabled.
  %&-line parsing enabled.
@@ -1032,7 +1032,7 @@ s/cm/cmsl10.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/c
 msy10.pfb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmsy7.p
 fb></usr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmti10.pfb></u
 sr/share/texlive/texmf-dist/fonts/type1/public/amsfonts/cm/cmtt10.pfb>
-Output written on thesis.pdf (8 pages, 226902 bytes).
+Output written on thesis.pdf (8 pages, 225772 bytes).
 PDF statistics:
  197 PDF objects out of 1000 (max. 8388607)
  168 compressed objects within 2 object streams

thesis/thesis.tex

@@ -295,7 +295,7 @@
 	
 	\section{Base Network}
 	\label{appendix:base}
-	Throughout this paper we rely on existing networks obtained via tomographic scans to serve as base for \textbf{generation} of larger networks and \textbf{simulation} of network flow. Their statistics are detailed in the figures below.
+	Throughout this paper we rely on existing networks obtained via tomographic scans to serve as a basis for \textbf{generation} of larger networks and \textbf{simulation} of network flow. Their statistics are detailed in the figures below.
 	\newpage
 	
 	\centering

visualize/print_pores.py

 import numpy as np
 
-nCells = [6, 6]
-
 def cell_pos(pos):
     from math import floor
     return [floor(pos[i] + 3.) for i in range(2)]