To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit c5a13ba7 authored by sfritschi's avatar sfritschi
Browse files

Added comparison function

parent 39ec32db
import sys
sys.path.append('../')
import matplotlib.pyplot as plt
import netflow.netgen as netgen_par
import netflow_serial.netgen as netgen_ser
from time import perf_counter
def compute_target_sizes(max_power2 = 2):
power = 0
d = 3
target_sizes = [[1, 1, 1]]
while (power < max_power2):
prev = target_sizes[d * power]
for i in range(d):
t = prev.copy()
t[i] <<= 1
prev = t
target_sizes.append(t)
power += 1
return target_sizes
def main():
if (len(sys.argv) != 2):
raise AssertionError("Usage: python3 compare.py <nthreads>")
basenet = netgen_par.load_network_from('../network/network.h5')
print("Network statistics:")
print("Throats: %d" % len(basenet.throats))
pore_throat_counts = [len(pore.throats) for pore in basenet.pores]
print("Max. number of throats per pore: %d" % max(pore_throat_counts))
print("Min. number of throats per pore: %d" % min(pore_throat_counts))
print("Avg. number of throats per pore: %f" % (sum(pore_throat_counts) / len(pore_throat_counts)))
nthreads = int(sys.argv[1])
print("Using: {} threads".format(nthreads))
targetsizes = compute_target_sizes(1)
cutoff = 0.5 * max([basenet.ub[i] - basenet.lb[i] \
for i in range(len(basenet.ub))])
speedups = []
sizes = []
times_par = []
times_ser = []
for i, target in enumerate(targetsizes):
sizes.append(netgen_par.prod(target))
print("Multiplicative generated size: {}".format(sizes[i]))
start = perf_counter()
dendro = netgen_par.generate_dendrogram(basenet=basenet, targetsize=target, \
cutoff=cutoff, sd=42, nthreads=nthreads, mute=True)
end = perf_counter()
elapsed = end - start
times_par.append(elapsed)
print("Elapsed par: %e" % elapsed)
print("Generated network statistics:")
print("Throats: %d" % len(dendro.throats))
print("Pores: %d" % len(dendro.pores))
start = perf_counter()
dendro = netgen_ser.generate_dendrogram(basenet=basenet, targetsize=target, \
cutoff=cutoff, sd=42, mute=True)
end = perf_counter()
elapsed = end - start
times_ser.append(elapsed)
speedups.append(times_ser[i] / times_par[i])
print("\nElapsed ser: %e" % elapsed)
print("Generated network statistics:")
print("Throats: %d" % len(dendro.throats))
print("Pores: %d" % len(dendro.pores))
print("")
plt.figure()
plt.title(r"Comparison Between Parallel and Serial Methods")
plt.xlabel(r"Multiplicative generated size")
plt.ylabel(r"Time [s]")
plt.plot(sizes, times_ser, '--rx', label="Serial")
plt.plot(sizes, times_par, '--bo', label="Parlallel using {} threads".format(nthreads))
plt.legend()
plt.grid(True)
plt.savefig('plots/comparison_time.png')
plt.close()
plt.figure()
plt.title("Speedup of Parallel Method over Serial Method")
plt.xlabel(r"Multiplicative generated size")
plt.ylabel(r"Speedup")
plt.plot(sizes, speedups, '--bo')
plt.grid(True)
plt.savefig('plots/comparison_speedup.png')
plt.close()
if __name__ == '__main__':
main()
......@@ -12,7 +12,7 @@ def main():
nthreads = int(sys.argv[1])
print("Using %d thread(s)" % nthreads)
basenet = netflow.load_network_from('../netflow/network/network.h5')
basenet = netflow.load_network_from('../network/network.h5')
print("Network statistics:")
pore_throat_counts = [len(pore.throats) for pore in basenet.pores]
print("Max. number of throats per pore: %d" % max(pore_throat_counts))
......
......@@ -10,7 +10,7 @@ def main():
if len(sys.argv) != 6:
raise AssertionError("Usage: python3 perf.py <2^min threads> <2^max threads> <n1> <n2> <n3>")
basenet = netflow.load_network_from('../netflow/network/network.h5')
basenet = netflow.load_network_from('../network/network.h5')
print("Network statistics:")
print("Throats: %d" % len(basenet.throats))
pore_throat_counts = [len(pore.throats) for pore in basenet.pores]
......@@ -54,7 +54,10 @@ def main():
plt.title(r"Strong Scaling for targetsize {}".format(target))
plt.xlabel(r"#processes")
plt.ylabel(r"Speedup")
plt.plot(threads, speedups, '--bo')
plt.plot(threads, speedups, '--bo', label=r"true speedup")
plt.plot(threads, threads, '-r', label=r"ideal speedup")
plt.grid(True)
plt.legend()
plt.savefig('plots/strong_scaling.png')
plt.close()
......
gen/plots/strong_scaling.png

21.2 KB | W: | H:

gen/plots/strong_scaling.png

25 KB | W: | H:

gen/plots/strong_scaling.png
gen/plots/strong_scaling.png
gen/plots/strong_scaling.png
gen/plots/strong_scaling.png
  • 2-up
  • Swipe
  • Onion skin
......@@ -163,8 +163,8 @@ class CellList:
def match_maker(self, pores: List[Pore], tid: int):
n = len(pores)
for i, pore in enumerate(pores):
if (tid == 0):
print(f"progress {((i+1) / n)*100.:.1f}%", end="\r", flush=True)
#if (tid == 0):
# print(f"progress {((i+1) / n)*100.:.1f}%", end="\r", flush=True)
candidates = self.find_candidates(pore)
# Find best candidate for each throat of pore
......@@ -950,12 +950,13 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
del ttbrs
# DEBUG
percent = nUnrealized / totalThroats * 100.
print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.3f}%)".\
format(len(throats), nUnrealized, percent))
avg_throat_diff /= count
print("Avg. throat length difference: %e" % avg_throat_diff)
print("Relative to Lmax: %.1f%%" % (avg_throat_diff / basenet.Lmax * 100.))
if (not mute):
percent = nUnrealized / totalThroats * 100.
print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.3f}%)".\
format(len(throats), nUnrealized, percent))
avg_throat_diff /= count
print("Avg. throat length difference: %e" % avg_throat_diff)
print("Relative to Lmax: %.1f%%" % (avg_throat_diff / basenet.Lmax * 100.))
# assemble and return network
network = Network(lb=[0.0 for k in range(d)],
......
__all__ = ['netgen', 'netflow']
# flow solution of periodic or in-/outflow pore-networks
#
# Daniel W. Meyer
# Institute of Fluid Dynamics, ETH Zurich
# January 2019
from typing import List, Set, Dict, Tuple # for type hints in argument lists
from .netgen import *
def solve_flow_periodic(network: Network, mu: float, c: int, \
P2L: float, tol: float = 1e-20, maxiter: int = 1e3, resinfo: bool = False) \
-> Tuple[Dict[Pore,float],Dict[Throat,float]]:
"""Solve flow problem for periodic network.
A mean pressure drop P/L [Pa/m] is applied in the c-direction with c>0.
L is the domain length and mu is the dynamic viscosity in [Pa*s = kg/(m*s)].
The pressure and flow solution is computed in a connected part of the pore
network that contains at least half of the pores in the network.
If no such network exists, an error is triggered.
Fluctuating pore pressures [Pa] and throat fluxes [m^3/s] are returned.
The first return value is a pore/pressure dictionary and the second a
throat/flux dictionary. Fluxes going out of pore1 of a throat are positive.
"""
L = [ub - lb for ub, lb in zip(network.ub,network.lb)]
# initialize
iterator = iter(network.pores)
seed = next(iterator)
tree = __grow({seed})
# find sufficiently large tree of connected pores
while len(tree) < len(network.pores)/2:
# find seed outside of tree being too small
for seed in iterator:
if (seed not in tree): break
else: raise AssertionError('Fragmented network!')
# grow tree based on new seed
tree = __grow({seed})
# assemble sparse pressure matrix
pore_dict = {pore:k for k, pore in enumerate(tree)}
A = {}; b = {} # pressure matrix and right hand side vector
for p_i, i in pore_dict.items(): # loop over pores i
if i == 0: # boundary condition
A[(0,0)] = 1.0; b[0] = 0.0
continue
for t_ij in p_i.throats:
# pressure system matrix A
T,p_j_pos = __add_coefficient_to(A, t_ij, pore_dict, p_i, i, mu, L)
# right hand side b
__add_to(b, i, T * (p_i.pos[c-1] - p_j_pos[c-1]) * P2L)
p = __amg_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
# calculate resulting fluxes
pm = lambda x: P2L*L[c-1] - (x[c-1] - network.lb[c-1]) * P2L # press. drop
Q = __calculate_fluxes(network, pore_dict, p, mu, L, pm)
return p, Q
def solve_flow_inout(network: Network, pin: float, pout: float, \
inpores: Set[Pore], outpores: Set[Pore], \
mu: float, tol: float = 1e-20, maxiter: int = 1e4, resinfo: bool = False) \
-> Tuple[Dict[Pore,float],Dict[Throat,float]]:
"""Solve flow problem for in-/outflow netw. with possibly periodic throats.
pin and pout [Pa] are applied at the in- and out-pores.
mu is the dynamic viscosity in [Pa*s = kg/(m*s)].
In case the sets of in-/out-pores are not connected, an error is triggered.
Pore pressures [Pa = kg/(m*s^2)] and throat fluxes [m^3/s] are returned.
The first return value is a pore/pressure dictionary and the second a
throat/flux dictionary. Fluxes going out of pore1 of a throat are positive.
"""
L = [ub - lb for ub, lb in zip(network.ub,network.lb)]
# are in- and outflow pores connected?
tree = __grow(inpores)
if len(tree & outpores) == 0:
raise AssertionError('In-/outflow disconnected!')
tree = tree | __grow(outpores) # collect pores connected to in-\outflow
# assemble sparse pressure matrix
pore_dict = {pore:k for k, pore in enumerate(tree)}
A = {}; b = {} # pressure matrix and right hand side vector
for p_i, i in pore_dict.items(): # loop over pores i
if p_i in inpores: # inflow boundary condition
A[(i,i)] = 1.0; b[i] = pin
continue
elif p_i in outpores: # outflow boundary condition
A[(i,i)] = 1.0; b[i] = pout
continue
for t_ij in p_i.throats: # pressure matrix coefficients
T,p_j_pos = __add_coefficient_to(A, t_ij, pore_dict, p_i, i, mu, L)
try:
p = __amg_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
#p = __cg_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
except:
print("Error in AMG solver: switching to ILU-GMRES solver")
p = __ilu_gmres_pressure_system(pore_dict, A, b, tol, maxiter, resinfo)
# calculate resulting fluxes
pm = lambda x: 0.0 # no mean pressure drop
Q = __calculate_fluxes(network, pore_dict, p, mu, L, pm)
return p, Q
def __amg_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
tol: float, maxiter: int, resinfo: bool = False) -> Dict[Pore,float]:
"""Solve pressure system with iterative algebraic multi-grid solver."""
import numpy as np
from scipy.sparse import csr_matrix
n = len(pore_dict)
# A of Ax = b
data = []; row = []; col = []
for ij in A:
data.append(A[ij]); row.append(ij[0]); col.append(ij[1])
An = csr_matrix((data, (row, col)), shape=(n, n))
# b
data = np.zeros((n,1))
for j in b: data[j] = b[j]
bn = data
# iterative solver (tol is relative residual r[k]/r[0] tolerance)
import pyamg
ml = pyamg.ruge_stuben_solver(An) # construct multigrid hierarchy
if resinfo:
rl = list(); j = 0 # residual list, iteration counter
def printres(xk): # residual info callback function
nonlocal rl, j; j = j+1
print(str(j) + ' ' + str(rl[-1]))
p = ml.solve(bn.T, tol=tol, maxiter=maxiter, cycle='W', \
callback=printres, residuals=rl)
else:
p = ml.solve(bn.T, tol=tol, maxiter=maxiter, cycle='W')
return {pore:p[k] for pore, k in pore_dict.items()}
def __cg_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
tol: float, maxiter: int, resinfo: bool = False) -> Dict[Pore,float]:
"""Solve pressure system with conjugate gradient solver."""
import numpy as np
from scipy.sparse import csc_matrix
n = len(pore_dict)
# A of Ax = b
data = []; row = []; col = []
for ij in A:
data.append(A[ij]); row.append(ij[0]); col.append(ij[1])
An = csc_matrix((data, (row, col)), shape=(n, n))
# b
data = np.zeros((n,1))
for j in b: data[j] = b[j]
bn = data
# iterative solver
from scipy.sparse import linalg
if resinfo:
j = 0 # iteration counter
def printres(xk): # residual info callback function
nonlocal An, bn, j; j = j+1
res = np.linalg.norm(An.dot(xk) - np.transpose(bn))
print(str(j) + ' ' + str(res))
p = linalg.cg(An, bn, tol=tol, maxiter=maxiter, callback=printres)[0]
else:
p = linalg.cg(An, bn, tol=tol, maxiter=maxiter)[0]
return {pore:p[k] for pore, k in pore_dict.items()}
def __ilu_gmres_pressure_system(pore_dict: Dict[Pore,int], A: Dict, b: Dict, \
tol: float, maxiter: int, resinfo: bool = False) -> Dict[Pore,float]:
"""Solve pressure system with incomplete LU precondit. + GMRES solver."""
import numpy as np
from scipy.sparse import csc_matrix
n = len(pore_dict)
# A of Ax = b
data = []; row = []; col = []
for ij in A:
data.append(A[ij]); row.append(ij[0]); col.append(ij[1])
An = csc_matrix((data, (row, col)), shape=(n, n))
# b
data = np.zeros((n,1))
for j in b: data[j] = b[j]
bn = data
# incomplete LU preconditioner
from scipy.sparse import linalg
invA_approx = linalg.spilu(An)
M = linalg.LinearOperator((n,n), invA_approx.solve)
# iterative solver
if resinfo:
j = 0 # iteration counter
def printres(xk): # residual info callback function
nonlocal j; j = j+1
import inspect
frame = inspect.currentframe().f_back
print(str(j) + ' ' + str(frame.f_locals['resid']))
p = linalg.gmres(An, bn, M=M, tol=tol, maxiter=maxiter, \
callback=printres)[0]
else:
p = linalg.gmres(An, bn, M=M, tol=tol, maxiter=maxiter)[0]
return {pore:p[k] for pore, k in pore_dict.items()}
def __add_to(D: Dict, i, d: float):
"""Append key i with value d to dict D or add v to D[i] if key i exists."""
if i in D: D[i] = D[i] + d # add to existing entry
else: D[i] = d # store new entry
def __add_coefficient_to(A: Dict, t_ij: Throat, pore_dict: Dict, \
p_i: Pore, i: int, mu: float, L: List[float]) -> Tuple[float,float]:
"""Add or modify coefficients in matrix of pressure system."""
# identify second pore j
p_j = t_ij.pore1
if p_i == p_j: p_j = t_ij.pore2
j = pore_dict[p_j] # index of pore j
# determine position of possibly periodic pore j
p_j_pos = p_j.pos.copy() # copy position
if len(t_ij.label) != 0: # is pore j a periodic pore?
direct = [int(k) for k in t_ij.label.split()] # periodicity
p_j_pos = [p_j.pos[k] + direct[k] * \
(+1 if t_ij.pore1 == p_i else -1) * L[k] \
for k in range(len(p_j.pos))] # shifted position
# determine and set coefficients in pressure matrix and r.h.s.
Lt = distance(p_i.pos, p_j_pos)
T = tc(t_ij.r, mu, Lt)
__add_to(A, (i,j), -T)
__add_to(A, (i,i), T)
return T, p_j_pos
def tc(r: float, mu: float, Lt:float):
"""Calculates the throat conductivity in [m^4*s/kg]."""
from math import pi
return pi/(4*32) * (2*r)**4 / (mu*Lt)
def __calculate_fluxes(network: Network, pore_dict: Dict[Pore,int], \
p: Dict[Pore,float], mu: float, L: List[float], pm) -> Dict[Throat,float]:
"""Calculate throat fluxes from pressure solution."""
Q = {}
for throat in network.throats:
p1 = throat.pore1; p2 = throat.pore2
if p1 not in pore_dict: continue # ignore disconnected pores
# pore positions needed to determine mean pressure drop
x1, x2 = throat_ends(throat, L)
# calculate flux
Lt = distance(x1, x2)
T = tc(throat.r, mu, Lt)
Q[throat] = T * ((p[p1]+pm(x1)) - (p[p2]+pm(x2)))
return Q
def __grow(seeds: Set[Pore]) -> Set[Pore]:
"""Starting from a seed set of pores, find all connected pores."""
# initialize
offsprings = seeds.copy(); tree = seeds.copy()
# grow tree
while len(offsprings) > 0: # loop over growth periods
branches = set() # initialize for next growth period
for node in offsprings.copy(): # loop over ...
for throat in node.throats: # ... next possible branches
branch = throat.pore1
if branch == node: branch = throat.pore2
if branch in tree: continue
# add next branch
branches.add(branch) # ... for new growth period
tree.add(branch) # ... to tree
offsprings = branches.copy()
return tree
def flux_balance(p: Dict[Pore,float], Q: Dict[Throat,float]) \
-> Dict[Pore,List[float]]:
"""Calculate the flux balance for each pore in the network.
The return value is a dictionary, where the keys are network pores and
where each value is a list of sums of outgoing, incoming, and all fluxes.
The fourth and last list element is the number of throat connections.
"""
Qb = {} # flux balance
# loop over connected pores
for pore in p.keys():
# add fluxes of connecting throats
Qp = 0; Qm = 0; Qsum = 0
for throat in pore.throats:
Qt = (+1 if pore == throat.pore1 else -1) * Q[throat]
if Qt > 0: Qp = Qp + Qt # outgoing fluxes
else: Qm = Qm + Qt # incoming fluxes
Qsum = Qsum + Qt # flux total
Qb[pore] = (Qp, Qm, Qsum, len(pore.throats))
return Qb
def flux_plane(network: Network, Q: Dict[Throat,float], \
c: int, x: float = None) -> float:
"""Total flux [m^3/s] in direction c over a balance plane at position x.
The balance plane position lb < x < ub and c > 0.
If x = None, x = (lb+ub)/2 is used.
"""
from math import copysign
Lc = network.ub[c-1]-network.lb[c-1]
# balance plane position
if x == None: x = (network.ub[c-1]+network.lb[c-1])/2
# sum up total flux
Qtot = 0.0
for throat, Qi in Q.items():
# throat end positions
x1 = throat.pore1.pos[c-1]
x2 = throat.pore2.pos[c-1]
# account for periodic throats
if len(throat.label) != 0:
direct = [int(k) for k in throat.label.split()]
x2 = x2 + direct[c-1]*Lc
if (x1 < x) and (x2 < x):
x1 = x1 + Lc; x2 = x2 + Lc
elif (x1 > x) and (x2 > x):
x1 = x1 - Lc; x2 = x2 - Lc
# add up flux of throats that are cut by balance plane
if (x2 - x) * (x - x1) > 0:
Qtot = Qtot + copysign(1,x2-x1)*Qi
return Qtot
def save_flow_to(filename: str, network: Network, \
p: Dict[Pore,float], Q: Dict[Throat,float]):
"""Save pressures and fluxes to hdf5 file."""
import h5py
f = h5py.File(filename, 'x')
# write network label
label = network.label.encode('ascii','ignore')
f.create_dataset('network_label', (1,), dtype='S'+str(len(label)),
data=[label])
# write pressure for each pore
p_grp = f.create_group('p')
pi = [0 for k in range(len(p))]; pp = [0.0 for k in range(len(p))]
for k, pore in enumerate(p):
pi[k] = pore.id; pp[k] = p[pore]
p_grp.create_dataset('pore', data=pi)
p_grp.create_dataset('p', data=pp)
# write flux for each throat
Q_grp = f.create_group('Q')
Qi = [0 for k in range(len(Q))]; QQ = [0.0 for k in range(len(Q))]
for k, throat in enumerate(Q):
Qi[k] = throat.id; QQ[k] = Q[throat]
Q_grp.create_dataset('throat', data=Qi)
Q_grp.create_dataset('Q', data=QQ)
def load_flow_from(filename: str, network: Network) \
-> Tuple[Dict[Pore,float],Dict[Throat,float]]:
"""Load pressures and fluxes from hdf5 file and link them with network.
To link the pore/pressure- and throat/flux-dictionaries with the
corresponding network objects, pore and throat ids are used.
"""
# translation of ids to pores and throats
pores = {pore.id:pore for pore in network.pores}
throats = {throat.id:throat for throat in network.throats}
# read file
import h5py
import numpy as np
f = h5py.File(filename, 'r')
# network label
label = f['network_label'][0].decode('utf-8')
if label != network.label:
raise AssertionError('Inconsistent network label!')
# read pressure for each pore
pi = np.array(f['p/pore']); pp = np.array(f['p/p'])
p = {}
for id, press in zip(pi,pp): p[pores[id]] = press
# read flux for each throat
Qi = np.array(f['Q/throat']); QQ = np.array(f['Q/Q'])
Q = {}
for id, flux in zip(Qi,QQ): Q[throats[id]] = flux
return p, Q
This diff is collapsed.
......@@ -26,3 +26,4 @@ cleanup:
rm -f *.so
rm -f cwrapper/wrapper.c
rm -rf ../netflow/__pycache__
rm -rf ../netflow_serial/__pycache__
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment