Commit 683437c3 authored by sfritschi's avatar sfritschi
Browse files

Consider entire neighborhood; Connector checks for periodic pores

parent aadd8dcf
...@@ -55,6 +55,7 @@ def main(): ...@@ -55,6 +55,7 @@ def main():
plt.ylabel(r"Speedup") plt.ylabel(r"Speedup")
plt.plot(threads, speedups, '--bo') plt.plot(threads, speedups, '--bo')
plt.savefig('plots/strong_scaling.png') plt.savefig('plots/strong_scaling.png')
plt.close()
if __name__ == '__main__': if __name__ == '__main__':
main() main()
...@@ -109,7 +109,7 @@ class CellList: ...@@ -109,7 +109,7 @@ class CellList:
""" """
def __init__(self, pores: List[Pore], copies: Dict[Pore, Set[Pore]], \ def __init__(self, pores: List[Pore], copies: Dict[Pore, Set[Pore]], \
domainSize: List[float], cellSize: float, \ domainSize: List[float], cellSize: float, \
basenet: Network, nInterior: int, isPeriodicCheck: bool): basenet: Network, nInterior: int):
from math import ceil from math import ceil
self.dim = len(domainSize) self.dim = len(domainSize)
self.Lmax = basenet.Lmax self.Lmax = basenet.Lmax
...@@ -131,53 +131,11 @@ class CellList: ...@@ -131,53 +131,11 @@ class CellList:
offset = 0 offset = 0
# Used to keep track of already realized throats # Used to keep track of already realized throats
self.poreThroatIndices = [set() for _ in range(nInterior)] self.poreThroatIndices = [set() for _ in range(nInterior)]
# Additional information needed for correctness if targetsize contains 1
self.isPeriodicCheck = isPeriodicCheck
if isPeriodicCheck:
# Maps periodic copy index to respective index inside domain
# {int: int}
self.originalIdx = {}
# Maps index from inside domain to List of all its periodic copies
# {int: List[int]}
self.periodicCopies = {}
L = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)] L = [ub-lb for lb,ub in zip(basenet.lb,basenet.ub)]
# Sort pores according to cell membership # Sort pores according to cell membership
self.poresSorted = [set() for _ in range(self.totalCells)] self.poresSorted = [set() for _ in range(self.totalCells)]
if isPeriodicCheck:
for poreIdx, pore in enumerate(pores):
cellIdx = self.pore_to_index(pore)
self.poresSorted[cellIdx].add(poreIdx)
self.poresPos += pore.pos
if poreIdx < nInterior:
# Pre-fix sum
nthroats = len(pore.throats)
offset += nthroats
self.throatOffsets.append(offset)
self.originalIdx[poreIdx] = poreIdx # identity
copyIndices = []
for cpore in copies[pore]:
copyIndices.append(cpore.index)
self.periodicCopies[poreIdx] = copyIndices
for throat in pore.throats:
self.throatL.append(distance(*throat_ends(throat, L)))
self.throatR.append(throat.r)
self.poreThroatIndices[poreIdx].add(throatIdx)
throatIdx += 1
else: # periodic copy of pore in interior
# Determine index of original pore
originalIdx, _ = pore.label.split(' ',1)
self.originalIdx[poreIdx] = int(originalIdx)
else:
for poreIdx, pore in enumerate(pores): for poreIdx, pore in enumerate(pores):
cellIdx = self.pore_to_index(pore) cellIdx = self.pore_to_index(pore)
...@@ -224,32 +182,20 @@ class CellList: ...@@ -224,32 +182,20 @@ class CellList:
def match_maker(self, pores: List[int], tid: int): def match_maker(self, pores: List[int], tid: int):
n = len(pores) n = len(pores)
for poreIdx in pores: for poreIdx in pores:
if (tid == 0):
print(f"progress {((poreIdx+1) / n)*100.:.1f}%", end="\r", flush=True)
candidates = self.find_candidates(poreIdx) candidates = self.find_candidates(poreIdx)
throats = self.fetch_throat_indices(poreIdx) throats = self.fetch_throat_indices(poreIdx)
for idx, throatIdx in enumerate(throats): for idx, throatIdx in enumerate(throats):
lt = self.throatL[throatIdx] lt = self.throatL[throatIdx]
# Find best match among candidates # Find best match among candidates
matchIdx, _ = min(candidates.items(), key=lambda il: abs(il[1] - lt)) matchIdx, _ = min(candidates.items(), key=lambda il: abs(il[1] - lt))
# Remove match from candidates # Remove match from candidates
del candidates[matchIdx] del candidates[matchIdx]
# Put best match index in poreThroatTable # Put best match index in poreThroatTable
self.set_throat(poreIdx, matchIdx, idx) self.set_throat(poreIdx, matchIdx, idx)
# Remove periodic copies as well if targetsize contains 1
if self.isPeriodicCheck:
# Get index of original pore in interior of domain
originalIdx = self.originalIdx[matchIdx]
if originalIdx in candidates:
del candidates[originalIdx]
for cporeIdx in self.periodicCopies[originalIdx]:
if cporeIdx in candidates:
del candidates[cporeIdx]
# Compute dictionary of all nbor candidates of given pore # Compute dictionary of all nbor candidates of given pore
def find_candidates(self, poreIdx: int) -> Dict[int, float]: def find_candidates(self, poreIdx: int) -> Dict[int, float]:
...@@ -294,16 +240,25 @@ class CellList: ...@@ -294,16 +240,25 @@ class CellList:
ub = self.dim + lb ub = self.dim + lb
return self.poresPos[lb:ub] return self.poresPos[lb:ub]
# Only consider subset of all neighboring cells (closest 7) # Consider all 27 neighboring cells
def nbor_indices(self, poreIdx: int) -> List[int]: def nbor_indices(self, poreIdx: int) -> List[int]:
cellIdx = self.pore_to_triplet(poreIdx) cellIdx = self.pore_to_triplet(poreIdx)
# i, j, k; i +- 1, j, k; i, j +- 1, k; i, j, k +- 1;
base = self.flatten(*cellIdx) base = self.flatten(*cellIdx)
widthx = self.nCells[0] widthx = self.nCells[0]
widthy = widthx * self.nCells[1] widthy = widthx * self.nCells[1]
return [base, base - 1, base + 1, \ return [base, base - 1, base + 1,
base - widthx, base + widthx, \ base - widthx, base + widthx,
base - widthy, base + widthy] base - widthy, base + widthy,
base + widthx - 1, base + widthx + 1,
base - widthx - 1, base - widthx + 1,
base + widthy - 1, base + widthy + 1,
base - widthy - 1, base - widthy + 1,
base + widthx + widthy, base + widthx + widthy - 1,
base + widthx + widthy + 1, base + widthx - widthy,
base + widthx - widthy - 1, base + widthx - widthy + 1,
base - widthx + widthy, base - widthx + widthy - 1,
base - widthx + widthy + 1, base - widthx - widthy,
base - widthx - widthy - 1, base - widthx - widthy + 1]
def nbor_cell_counts(self, neighborhood: List[int]) -> List[int]: def nbor_cell_counts(self, neighborhood: List[int]) -> List[int]:
return [len(self.poresSorted[i]) for i in neighborhood] return [len(self.poresSorted[i]) for i in neighborhood]
...@@ -791,17 +746,14 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \ ...@@ -791,17 +746,14 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
n = len(pores) n = len(pores)
copies = __add_buffer_layers(pores, L, basenet.Lmax) copies = __add_buffer_layers(pores, L, basenet.Lmax)
# Flag for handling the case where original nbor pore and its periodic
# copies are potential candidates of a given pore at the same time
# (necessary only if target size is 1 in at least one dimension)
isPeriodicCheck = 1 in targetsize isPeriodicCheck = 1 in targetsize
# Domain size including periodic buffer layers on all sides # Domain size including periodic buffer layers on all sides
trueDomainSize = [L[i] + 2 * basenet.Lmax for i in range(d)] trueDomainSize = [L[i] + 2 * basenet.Lmax for i in range(d)]
# Initialize cell-list; Place each pore in respective cell-set # Initialize cell-list; Place each pore in respective cell-set
cellList = CellList(pores, copies, trueDomainSize, basenet.Lmax, \ cellList = CellList(pores, copies, trueDomainSize, basenet.Lmax, \
basenet, n, isPeriodicCheck) basenet, n)
# Evenly distribute pores among threads # Evenly distribute pores (in interior) among threads
remainder = n % nthreads remainder = n % nthreads
loads = [n // nthreads] * nthreads loads = [n // nthreads] * nthreads
for i in range(remainder): for i in range(remainder):
...@@ -845,6 +797,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \ ...@@ -845,6 +797,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
matches = cellList.fetch_matches(poreIdx) matches = cellList.fetch_matches(poreIdx)
throatIndices = cellList.fetch_throat_indices(poreIdx) throatIndices = cellList.fetch_throat_indices(poreIdx)
conporesInd = set() # Set of connected pores (needed if target size contains 1)
for i, match in enumerate(matches): for i, match in enumerate(matches):
throatIdx = throatIndices[i] throatIdx = throatIndices[i]
...@@ -864,10 +817,12 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \ ...@@ -864,10 +817,12 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
nborTtbr = cellList.poreThroatIndices[nbor.index] nborTtbr = cellList.poreThroatIndices[nbor.index]
total += 1 total += 1
if (len(nborTtbr) == 0): if (len(nborTtbr) == 0 or nbor.index in conporesInd):
n_already_taken += 1 n_already_taken += 1
continue continue
conporesInd.add(nbor.index)
r = cellList.throatR[throatIdx] r = cellList.throatR[throatIdx]
targetLen = cellList.throatL[throatIdx] targetLen = cellList.throatL[throatIdx]
# Compute actual throat length between pores # Compute actual throat length between pores
...@@ -993,7 +948,6 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \ ...@@ -993,7 +948,6 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
del cellList del cellList
# DEBUG # DEBUG
"""
percent = k / (k + len(throats)) * 100. percent = k / (k + len(throats)) * 100.
print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.1f}%)".\ print("\b"*24 + "{0:d} throats in total, {1:d} unrealised ({2:.1f}%)".\
format(len(throats), k, percent)) format(len(throats), k, percent))
...@@ -1001,7 +955,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \ ...@@ -1001,7 +955,7 @@ def generate_dendrogram(basenet: Network, targetsize: List[int], \
print("Avg. throat length difference: %e" % avg_throat_diff) print("Avg. throat length difference: %e" % avg_throat_diff)
print("Relative to Lmax: %.1f%%" % (avg_throat_diff / basenet.Lmax * 100.)) print("Relative to Lmax: %.1f%%" % (avg_throat_diff / basenet.Lmax * 100.))
print("Percentage where match was fully-connected: %.1f%%" % (n_already_taken / total * 100.)) print("Percentage where match was fully-connected: %.1f%%" % (n_already_taken / total * 100.))
"""
# assemble and return network # assemble and return network
network = Network(lb=[0.0 for k in range(d)], network = Network(lb=[0.0 for k in range(d)],
ub=L, Lmax=basenet.Lmax, label='from_' + basenet.label) ub=L, Lmax=basenet.Lmax, label='from_' + basenet.label)
......
Supports Markdown
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