'''
Split blocks is combination of help from Dave Rigby and Tim Beach
'''
from .face import Face
from .block import Block
from typing import List
from enum import Enum
from math import gcd, sqrt
import numpy as np
[docs]
class Direction(Enum):
i = 0
j = 1
k = 2
[docs]
def max_aspect_ratio(X:np.ndarray,Y:np.ndarray,Z:np.ndarray,ix:int,jx:int,kx:int):
"""Finds the maximum cell aspect ratio
Args:
X (np.ndarray): 3 dimensional Array containing X values. X[i,j,k]
Y (np.ndarray): 3 dimensional Array containing X values. X[i,j,k]
Z (np.ndarray): 3 dimensional Array containing X values. X[i,j,k]
Returns:
float: Maximum value of aspect ratio
"""
[ix,jx,kx] = X.shape # ix, jx, kx are the max values like IMAX, JMAX, KMAX
ix-=1; jx-=1; kx-=1 # Python is 0 index so need to subtract 1 from ix, iy, kx to reference the last value
# Each column is a corner
# [corner 1, corner 2, corner 3, corner 4]
i1 = [0, ix, 0, 0]
j1 = [0, 0, jx, 0]
k1 = [0, 0, 0, kx]
# Appending more corners []
i1 = i1 + [0,ix,ix,ix]
j1 = j1 + [jx,0,jx,jx]
k1 = k1 + [kx,kx,0,kx]
# i2, j2, k2 are the first nodes diagonally from the block corners
i2 = [1, ix-1, 1, 1]
j2 = [1, 1, jx-1, 1]
k2 = [1, 1, 1, kx-1]
i2 = i2 + [1, ix-1, ix-1, ix-1]
j2 = j2 + [jx-1, 1, jx-1, jx-1]
k2 = k2 + [kx-1, kx-1, 1, kx-1]
ds = list()
for n in range(len(i2)):
ds.append(sqrt((X(i2[n],j1[n],k1[n])-X(i1[n],j1[n],k1[n]))**2 +
(Y(i2[n],j1[n],k1[n])-Y(i1[n],j1[n],k1[n]))**2 +
(Z(i2[n],j1[n],k1[n])-Z(i1[n],j1[n],k1[n]))**2)
)
ds.append(sqrt((X(i1[n],j2[n],k1[n])-X(i1[n],j1[n],k1[n]))**2 +
(Y(i1[n],j2[n],k1[n])-Y(i1[n],j1[n],k1[n]))**2 +
(Z(i1[n],j2[n],k1[n])-Z(i1[n],j1[n],k1[n]))**2)
)
ds.append(sqrt((X(i1[n],j1[n],k2[n])-X(i1[n],j1[n],k1[n]))**2 +
(Y(i1[n],j1[n],k2[n])-Y(i1[n],j1[n],k1[n]))**2 +
(Z(i1[n],j1[n],k2[n])-Z(i1[n],j1[n],k1[n]))**2)
)
aspect = [0,0,0]
if ds[0]>0:
aspect[0] = max(ds[1],ds[2])/ds[0]
elif(ds[1]>0):
aspect[1] = max(ds[2],ds[0])/ds[1]
elif(ds[2]>0):
aspect[2] = max(ds[0],ds[1])/ds[2]
return max(aspect)
def __step_search(total_cells:int,greatest_common_divisor:int,ncells_per_block:int,denominator:float,direction:str='forward'):
"""Searches for the right step size that leads to block splits with the same greatest common divisor. Same greatest common denominator means you will always have the same multi-grid when you solve.
Args:
total_cells (int): [description]
greatest_common_divisor (int): [description]
ncells_per_block (int): [description]
denominator (float): If the direction is "i" then denominator is JMAX*KMAX
direction (str, optional): 'forward' or 'backward' This chooses to increment step direction by +1 or -1. Defaults to 'forward'.
Returns:
[type]: [description]
"""
step_size = round(ncells_per_block/denominator) # initial Guess
initial_guess = step_size
Number_of_Cells_Remaining = total_cells % (step_size*denominator)
IJK_MAX_Remainder = Number_of_Cells_Remaining/denominator
increment = -1
if direction.lower() == 'forward':
increment=1
while ( (step_size %greatest_common_divisor != 0) or
((IJK_MAX_Remainder-1)%greatest_common_divisor !=0) and
step_size>initial_guess/2 and
step_size<initial_guess*1.5):
if (step_size %greatest_common_divisor == 0) and (IJK_MAX_Remainder-1)%greatest_common_divisor ==0:
break
step_size+=increment
Number_of_Cells_Remaining = total_cells % ((step_size)*denominator)
IJK_MAX_Remainder = Number_of_Cells_Remaining/denominator
if step_size %greatest_common_divisor != 0: # This checks if each of the split blocks is divisible by gcd
if (IJK_MAX_Remainder-1)%greatest_common_divisor !=0: # This checks the remaining/final block to see if it is divisible by gcd
return -1
return step_size
[docs]
def split_blocks(blocks:List[Block], ncells_per_block:int,direction:Direction=None):
"""Split blocks is used to divide an array of blocks based on number of cells per block. This code maintains the greatest common denominator of the parent block. Number of cells per block is simply an estimate of how many you want. The actual number will change to meet the greatest common denominator (GCD). GCD of 4 means multigrid of 3 e.g. grid/4 (coarse), 2 (fine), and 1 (finest). If a direction is not specified then for each block the longest index either i,j, or k is used.
Wisdom from Dave Rigby:
For example, for radial equilibrium we must integrate across the span. Some codes (GlennHT used to) would want a single block across the entire span. In that case you would want some additional control.
Another example might be if you would like a block to include the entire boundary layer. In that case you might introduce an aspect ratio control.
Args:
blocks (List[Block]): List of blocks
ncells_per_block (int): number of cells desired per block
direction (Direction): direction to split the blocks in. Direction.(i,j,k). Defaults to None. None means it will pick the direction for you based on which is greater IMAX, JMAX, or KMAX
Returns:
Blocks (List[Block]): list of blocks split in the specified direction
"""
direction_to_use = direction # store the user input variable
new_blocks = list()
for block_indx in range(len(blocks)):
block = blocks[block_indx]
total_cells = block.IMAX*block.JMAX*block.KMAX
if direction==None:
indx = np.argmin(np.array([block.IMAX,block.JMAX,block.KMAX]))
if indx == 0:
direction_to_use=Direction.i
elif indx == 1:
direction_to_use=Direction.j
elif indx == 2:
direction_to_use=Direction.k
if total_cells>ncells_per_block:
# Use greatest common divsor to maintain multi-grid so say the entire block is divisible by 4 then we want to maintain than for all the splits!
greatest_common_divisor =gcd(block.IMAX-1, gcd(block.JMAX-1, block.KMAX-1)) # Gets the maximum number of partitions that we can make for this given block
if direction_to_use == Direction.i:
# In order to get close to the number of cells per block, we need to control how many steps of the greatest_common_divisor to advance so for example if you have a multigrid mesh that has gcd of 16 (fine) => 8 (coarse) => 4 (coarser) => 2 (coarsest) and you want 400K cells per block then JMAX*KMAX*gcd*some_factor has to be close to 400K cells
denominator = block.JMAX*block.KMAX
step_size = __step_search(total_cells,greatest_common_divisor,ncells_per_block,denominator,direction='backward')
if step_size==-1:
step_size = __step_search(total_cells,greatest_common_divisor,ncells_per_block,denominator,direction='forward')
if step_size==-1:
assert('no valid step size found, do you have multi-block? gcd > 1')
# step_size-1 is the IMAX of the sub_blocks e.g. 0 to 92 this shows IMAX=93, (93-1) % 4 = 0 (good)
iprev = 0
for i in range(step_size,block.IMAX,step_size):
if (i+1) > block.IMAX:
break
X = block.X[iprev:i+1,:,:] # New X, Y, Z splits
Y = block.Y[iprev:i+1,:,:] # This indexes to iprev:i so if iprev=2 and i = 10 it will go from 2 to 9
Z = block.Z[iprev:i+1,:,:]
iprev=i # Blocks have to share the same face, Pick the previous face
new_blocks.append(Block(X,Y,Z))
# Check for remainder
if i+1 < block.IMAX:
# Add remainder to last block
X = block.X[i:,:,:] # New X, Y, Z splits
Y = block.Y[i:,:,:]
Z = block.Z[i:,:,:]
new_blocks.append(Block(X,Y,Z))
elif direction_to_use == Direction.j:
denominator = block.IMAX*block.KMAX
step_size = __step_search(total_cells,greatest_common_divisor,ncells_per_block,denominator,direction='backward')
if step_size==-1:
step_size = __step_search(total_cells,greatest_common_divisor,ncells_per_block,denominator,direction='forward')
if step_size==-1:
assert('no valid step size found, do you have multi-block? gcd > 1')
jprev = 0
for j in range(step_size,block.JMAX,step_size):
if (j+1) > block.IMAX:
break
X = block.X[:,jprev:j,:] # New X, Y, Z splits
Y = block.Y[:,jprev:j,:]
Z = block.Z[:,jprev:j,:]
jprev=j
new_blocks.append(Block(X,Y,Z))
# Check for remainder
if j+1 < block.JMAX:
# Add remainder to last block
X = block.X[:,j:,:] # New X, Y, Z splits
Y = block.Y[:,j:,:]
Z = block.Z[:,j:,:]
new_blocks.append(Block(X,Y,Z))
else:
denominator = block.IMAX*block.JMAX
step_size = __step_search(total_cells,greatest_common_divisor,ncells_per_block,denominator,direction='backward')
if step_size==-1:
step_size = __step_search(total_cells,greatest_common_divisor,ncells_per_block,denominator,direction='forward')
if step_size==-1:
assert('no valid step size found, do you have multi-block? gcd > 1')
kprev = 0
for k in range(step_size,block.KMAX,step_size):
if (k+1) > block.KMAX:
break
X = block.X[:,:,kprev:k+1] # New X, Y, Z splits
Y = block.Y[:,:,kprev:k+1]
Z = block.Z[:,:,kprev:k+1]
kprev=k
new_blocks.append(Block(X,Y,Z))
# Check for remainder
if k+1 < block.KMAX:
# Add remainder to last block
X = block.X[:,:,k:] # New X, Y, Z splits
Y = block.Y[:,:,k:]
Z = block.Z[:,:,k:]
new_blocks.append(Block(X,Y,Z)) # replace it
else:
new_blocks.append(block)
return new_blocks