from pumapy.utilities.workspace import Workspace
import math
import numpy as np
[docs]def size_check(size):
if not (isinstance(size, tuple) and len(size) == 3):
raise Exception("Error, invalid size, must be tuple of size 3")
if size[0] <= 0:
raise Exception("Error, invalid size, must be >= 0")
if size[1] <= 0:
raise Exception("Error, invalid size, must be >= 0")
if size[2] <= 0:
raise Exception("Error, invalid size, must be >= 0")
[docs]def range_inc(var, var_range, var_name):
if var < var_range[0] or var > var_range[1]:
raise Exception("Error, invalid " + str(var_name) + ", must be between " + str(var_range[0]) + " and " + str(
var_range[1]) + ", inclusive")
[docs]def range_exc(var, var_range, var_name):
if var <= var_range[0] or var >= var_range[1]:
raise Exception("Error, invalid " + str(var_name) + ", must be between " + str(var_range[0]) + " and " + str(
var_range[1]) + ", exclusive")
[docs]def greater_than_inc(var, lim, var_name):
if var < lim:
raise Exception("Error, invalid " + str(var_name) + ", must be greater than " + str(lim) + ", inclusive")
[docs]def greater_than_exc(var, lim, var_name):
if var <= lim:
raise Exception("Error, invalid " + str(var_name) + ", must be greater than " + str(lim) + ", exclusive")
[docs]def less_than_inc(var, lim, var_name):
if var > lim:
raise Exception("Error, invalid " + str(var_name) + ", must be less than " + str(lim) + ", inclusive")
[docs]def less_than_exc(var, lim, var_name):
if var >= lim:
raise Exception("Error, invalid " + str(var_name) + ", must be less than " + str(lim) + ", exclusive")
[docs]def check_ws_cutoff(workspace, cutoff):
if not isinstance(workspace, Workspace):
raise Exception("Need to pass pumapy.Workspace.")
if workspace.matrix.shape == [1, 1, 1]:
raise Exception("Empty workspace")
if len(cutoff) != 2:
raise Exception("Cutoff should be of length 2. Ex: (128,255)")
if min(cutoff[0], cutoff[1]) < 0:
raise Exception("Invalid cutoff. Must be positive")
if cutoff[0] > cutoff[1]:
raise Exception("Invalid cutoff. cutoff[0] should be <= cutoff[1]")
[docs]def estimate_max_memory(material_property, workspace_shape, solver_type='iterative', need_to_orient=False,
mf=True, perm_fluid_vf=1.):
""" Compute an estimate of the extra maximum memory required to run a specified material property on a domain
:param material_property: property to estimate, options:
'anisotropic_conductivity', 'isotropic_conductivity', 'tortuosity', 'elasticity', 'permeability'
:type material_property: string
:param workspace_shape: size of the domain to compute
:type workspace_shape: (int, int, int)
:param solver_type: type of solver, options: 'direct', 'iterative'
:type solver_type: string
:param need_to_orient: domain with orientation (needed for anisotropic conductivity and elasticity)
:type need_to_orient: bool
:param mf: using matrix-free construction of the system (needed for FE methods)
:type mf: bool
:return: number of Bytes
:rtype: int
"""
# missing properties: 'orientation', 'radiation'
mat_properties = ['anisotropic_conductivity', 'isotropic_conductivity', 'tortuosity', 'elasticity', 'elasticity_fe',
'permeability']
if material_property not in mat_properties:
raise Exception(f"material_property input can only be one of the following types: {mat_properties}")
def convert_bytes_size(size_bytes):
if size_bytes == 0:
return "0B"
size_name = ("B", "KB", "MB", "GB", "TB", "PB", "EB", "ZB", "YB")
i = int(math.floor(math.log(size_bytes, 1024)))
p = math.pow(1024, i)
s = round(size_bytes / p, 2)
return "%s %s" % (s, size_name[i])
# sizes in bytes
float_size = np.zeros(1, dtype=float).itemsize
uint32_size = np.zeros(1, dtype=np.uint32).itemsize
uint16_size = np.zeros(1, dtype=np.uint16).itemsize
len_x, len_y, len_z = workspace_shape
total_bytes = 0
# General idea: peak memory should be reached at Amat sparse matrix creation
if material_property in ["anisotropic_conductivity", "elasticity"]:
# because of the use of padded domains
len_x += 2
len_y += 2
len_z += 2
if material_property == "anisotropic_conductivity":
values_in_Amat_rows = 27
dof = 1. # degrees of freedom
else:
values_in_Amat_rows = 81
dof = 3.
len_xyz = len_x * len_y * len_z
# Amat size (V, I, J)
A_rows = dof * values_in_Amat_rows * len_xyz
total_bytes += 2 * (A_rows * float_size + 2 * A_rows * uint32_size)
# bvec size
bvec = dof * len_y * len_z * float_size + 2 * len_xyz * uint32_size
total_bytes += bvec
# ws_pad_size
total_bytes += len_xyz * uint16_size
# case dependent variables
if solver_type != "direct":
total_bytes += dof * len_xyz * float_size # initial guess size
if need_to_orient:
total_bytes += 3 * len_xyz * float_size # orient_pad_size
elif material_property in ["isotropic_conductivity", "tortuosity"]:
# Amat size (V, I, J)
A_rows = (len_x - 2) * len_y * len_z * 7 + 2 * len_y * len_z
total_bytes += 2 * (A_rows * float_size + 2 * A_rows * uint32_size)
len_xyz = len_x * len_y * len_z
# bvec size
total_bytes += len_xyz * float_size
# ws_pad_size
total_bytes += len_xyz * uint16_size
if solver_type != "direct":
total_bytes += len_xyz * float_size # initial guess size
elif material_property == "permeability":
n_elems = len_x * len_y * len_z
n_dofs = n_elems * 4 * perm_fluid_vf
n_nodes = (len_x + 1) * (len_y + 1) * (len_z + 1)
if mf == False or solver_type == 'direct':
total_bytes += (16 * n_elems + 64 * n_nodes + 2 * 64 * n_dofs + 18 * 64 * n_elems) / 8
else:
total_bytes += (16 * n_elems + 64 * n_nodes + 5 * 64 * n_dofs) / 8
elif material_property == "elasticity_fe":
n_elems = len_x * len_y * len_z
n_dofs = n_elems * 3
n_nodes = (len_x + 1) * (len_y + 1) * (len_z + 1)
if mf == False or solver_type == 'direct' or need_to_orient:
total_bytes += (16 * n_elems + 64 * n_nodes + 2 * 64 * n_dofs + 18 * 64 * n_elems) / 8
else:
total_bytes += (16 * n_elems + 64 * n_nodes + 5 * 64 * n_dofs) / 8
# elif material_property == "orientation":
# pass
#
# elif material_property == "radiation":
# pass
print(f"Approximate memory requirement for simulation: {convert_bytes_size(int(round(total_bytes)))}")
return total_bytes
[docs]def set_random_seed(seed):
""" Set random seed for scipy and numpy to make results reproducible.
NB: if you want to generate the same material twice in the same process, you need to call it twice (see example)
:param seed: random seed
:type seed: int
:Example:
>>> import pumapy as puma
>>> puma.set_random_seed(1)
>>> ws = puma.generate_random_spheres((100, 100, 100), 20, 0.5, allow_intersect=True, segmented=False)
Approximately ... spheres to be generated...
>>> # puma.render_volume(ws) # to visualize it
>>> puma.set_random_seed(1) # need to call it again to get the same domain!
>>> ws = puma.generate_random_spheres((100, 100, 100), 20, 0.5, allow_intersect=True, segmented=False)
Approximately ... spheres to be generated...
>>> # puma.render_volume(ws)
"""
np.random.seed(seed)