Skip to content

Commit

Permalink
make mesh.get_indices more compact
Browse files Browse the repository at this point in the history
  • Loading branch information
@ilhamv
ilhamv committed Aug 20, 2024
1 parent 62cca97 commit 19197fb
Show file tree
Hide file tree
Showing 5 changed files with 154 additions and 126 deletions.
4 changes: 1 addition & 3 deletions mcdc/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,9 +340,7 @@ def surface_distance(x, y, z, t, ux, uy, uz, v, surface, mcdc):
d_max = (t_max - t) * v

div = G * ux + H * uy + I_ * uz + J1 / v
distance = -evaluation / (
G * ux + H * uy + I_ * uz + J1 / v
)
distance = -evaluation / (G * ux + H * uy + I_ * uz + J1 / v)

# Go beyond current movement slice?
if distance > d_max:
Expand Down
28 changes: 14 additions & 14 deletions mcdc/iqmc/iqmc_kernel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -300,13 +300,13 @@ def iqmc_prepare_particles(mcdc):
P_new["ux"], P_new["uy"], P_new["uz"] = iqmc_sample_isotropic_direction(
lds[n, 1], lds[n, 5]
)
x = P_new['x']
y = P_new['y']
z = P_new['z']
t = P_new['t']
ux = P_new['ux']
uy = P_new['uy']
uz = P_new['uz']
x = P_new["x"]
y = P_new["y"]
z = P_new["z"]
t = P_new["t"]
ux = P_new["ux"]
uy = P_new["uy"]
uz = P_new["uz"]
ix, iy, iz, it, outside = mesh_.get_indices(x, y, z, t, ux, uy, uz, mesh)
q = Q[:, it, ix, iy, iz].copy()
dV = iqmc_cell_volume(ix, iy, iz, mesh)
Expand Down Expand Up @@ -347,13 +347,13 @@ def iqmc_score_tallies(P, distance, mcdc):
SigmaT = material["total"]
mat_id = P["material_ID"]

x_ = P['x']
y_ = P['y']
z_ = P['z']
t_ = P['t']
ux_ = P['ux']
uy_ = P['uy']
uz_ = P['uz']
x_ = P["x"]
y_ = P["y"]
z_ = P["z"]
t_ = P["t"]
ux_ = P["ux"]
uy_ = P["uy"]
uz_ = P["uz"]

x, y, z, t, outside = mesh_.get_indices(x_, y_, z_, t_, ux_, uy_, uz_, mesh)
if outside:
Expand Down
47 changes: 6 additions & 41 deletions mcdc/kernel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,16 +40,9 @@ def domain_crossing(P, mcdc):
if mcdc["technique"]["domain_decomposition"]:
mesh = mcdc["technique"]["dd_mesh"]
# Determine which dimension is crossed
x = P['x']
y = P['y']
z = P['z']
t = P['t']
ux = P['ux']
uy = P['uy']
uz = P['uz']
ix, iy, iz, it, outside = mesh_.get_indices(x, y, z, t, ux, uy, uz, mesh)

d_idx = mcdc['dd_idx']
ix, iy, iz, it, outside = mesh_.get_indices(P, mesh)

d_idx = mcdc["dd_idx"]
d_Nx = mcdc["technique"]["dd_mesh"]["x"].size - 1
d_Ny = mcdc["technique"]["dd_mesh"]["y"].size - 1
d_Nz = mcdc["technique"]["dd_mesh"]["z"].size - 1
Expand Down Expand Up @@ -423,15 +416,7 @@ def particle_in_domain(P, mcdc):
d_ix = int(d_idx - d_Nx * d_Ny * d_iz - d_Nx * d_iy)

mesh = mcdc["technique"]["dd_mesh"]
x = P['x']
y = P['y']
z = P['z']
t = P['t']
ux = P['ux']
uy = P['uy']
uz = P['uz']

x_cell, y_cell, z_cell, t_cell, outside = mesh_.get_indices(x, y, z, t, ux, uy, uz, mesh)
x_cell, y_cell, z_cell, t_cell, outside = mesh_.get_indices(P, mesh)

if d_ix == x_cell:
if d_iy == y_cell:
Expand Down Expand Up @@ -1685,17 +1670,8 @@ def score_tracklength(P, distance, mcdc):
tally = mcdc["tally"]
material = mcdc["materials"][P["material_ID"]]

# Particle coordinate
x = P["x"]
y = P["y"]
z = P["z"]
t = P["t"]
ux = P["ux"]
uy = P["uy"]
uz = P["uz"]

# Get indices
ix, iy, iz, it, outside = mesh_.get_indices(x, y, z, t, ux, uy, uz, tally["mesh"])
ix, iy, iz, it, outside = mesh_.get_indices(P, tally["mesh"])
mu, azi = mesh_get_angular_index(P, tally["mesh"])
g, outside_energy = mesh_get_energy_index(P, tally["mesh"], mcdc)

Expand Down Expand Up @@ -2816,19 +2792,8 @@ def branchless_collision(P, prog):
def weight_window(P, prog):
mcdc = adapt.device(prog)

# Particle coordinate
x = P["x"]
y = P["y"]
z = P["z"]
t = P["t"]
ux = P["ux"]
uy = P["uy"]
uz = P["uz"]

# Get indices
ix, iy, iz, it, outside = mesh_.get_indices(
x, y, z, t, ux, uy, uz, mcdc["technique"]["ww_mesh"]
)
ix, iy, iz, it, outside = mesh_.get_indices(P, mcdc["technique"]["ww_mesh"])

# Target weight
w_target = mcdc["technique"]["ww"][it, ix, iy, iz]
Expand Down
58 changes: 38 additions & 20 deletions mcdc/mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@


@njit
def get_indices(x, y, z, t, ux, uy, uz, mesh):
def get_indices(particle, mesh):
"""
Get mesh indices given the coordinate
Expand All @@ -15,29 +15,47 @@ def get_indices(x, y, z, t, ux, uy, uz, mesh):
# Check if coordinate is outside the mesh grid
outside = False
if (
x < mesh["x"][0] - COINCIDENCE_TOLERANCE
or x > mesh["x"][-1] + COINCIDENCE_TOLERANCE
or y < mesh["y"][0] - COINCIDENCE_TOLERANCE
or y > mesh["y"][-1] + COINCIDENCE_TOLERANCE
or z < mesh["z"][0] - COINCIDENCE_TOLERANCE
or z > mesh["z"][-1] + COINCIDENCE_TOLERANCE
or t < mesh["t"][0] - COINCIDENCE_TOLERANCE
or t > mesh["t"][-1] + COINCIDENCE_TOLERANCE
or (abs(x - mesh["x"][0]) < COINCIDENCE_TOLERANCE and ux < 0.0)
or (abs(x - mesh["x"][-1]) < COINCIDENCE_TOLERANCE and ux > 0.0)
or (abs(y - mesh["y"][0]) < COINCIDENCE_TOLERANCE and uy < 0.0)
or (abs(y - mesh["y"][-1]) < COINCIDENCE_TOLERANCE and uy > 0.0)
or (abs(z - mesh["z"][0]) < COINCIDENCE_TOLERANCE and uz < 0.0)
or (abs(z - mesh["z"][-1]) < COINCIDENCE_TOLERANCE and uz > 0.0)
or (abs(t - mesh["t"][-1]) < COINCIDENCE_TOLERANCE and 1.0 > 0.0)
particle["x"] < mesh["x"][0] - COINCIDENCE_TOLERANCE
or particle["x"] > mesh["x"][-1] + COINCIDENCE_TOLERANCE
or particle["y"] < mesh["y"][0] - COINCIDENCE_TOLERANCE
or particle["y"] > mesh["y"][-1] + COINCIDENCE_TOLERANCE
or particle["z"] < mesh["z"][0] - COINCIDENCE_TOLERANCE
or particle["z"] > mesh["z"][-1] + COINCIDENCE_TOLERANCE
or particle["t"] < mesh["t"][0] - COINCIDENCE_TOLERANCE
or particle["t"] > mesh["t"][-1] + COINCIDENCE_TOLERANCE
or (
abs(particle["x"] - mesh["x"][0]) < COINCIDENCE_TOLERANCE
and particle["ux"] < 0.0
)
or (
abs(particle["x"] - mesh["x"][-1]) < COINCIDENCE_TOLERANCE
and particle["ux"] > 0.0
)
or (
abs(particle["y"] - mesh["y"][0]) < COINCIDENCE_TOLERANCE
and particle["uy"] < 0.0
)
or (
abs(particle["y"] - mesh["y"][-1]) < COINCIDENCE_TOLERANCE
and particle["uy"] > 0.0
)
or (
abs(particle["z"] - mesh["z"][0]) < COINCIDENCE_TOLERANCE
and particle["uz"] < 0.0
)
or (
abs(particle["z"] - mesh["z"][-1]) < COINCIDENCE_TOLERANCE
and particle["uz"] > 0.0
)
or (abs(particle["t"] - mesh["t"][-1]) < COINCIDENCE_TOLERANCE and 1.0 > 0.0)
):
outside = True
return -1, -1, -1, -1, outside

ix = get_grid_index(x, ux, mesh["x"])
iy = get_grid_index(y, uy, mesh["y"])
iz = get_grid_index(z, uz, mesh["z"])
it = get_grid_index(t, 1.0, mesh["t"])
ix = get_grid_index(particle["x"], particle["ux"], mesh["x"])
iy = get_grid_index(particle["y"], particle["uy"], mesh["y"])
iz = get_grid_index(particle["z"], particle["uz"], mesh["z"])
it = get_grid_index(particle["t"], 1.0, mesh["t"])

return ix, iy, iz, it, outside

Expand Down
Loading

0 comments on commit 19197fb

Please sign in to comment.