Skip to content

Commit

Permalink
update precursor source generator
Browse files Browse the repository at this point in the history
  • Loading branch information
@ilhamv
ilhamv committed Nov 26, 2024
1 parent 5cc5117 commit 4d5293a
Show file tree
Hide file tree
Showing 2 changed files with 26 additions and 59 deletions.
82 changes: 24 additions & 58 deletions mcdc/loop.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,22 +340,11 @@ def loop_source(seed, data, mcdc):
work_end = work_start + work_size

for idx_work in range(work_size):
# =====================================================================
# Generate a source particle
# =====================================================================

generate_source_particle(work_start, idx_work, seed, mcdc)

# =====================================================================
# Run the source particle and its secondaries
# =====================================================================

exhaust_active_bank(data, mcdc)

# =====================================================================
# Closeout
# =====================================================================

source_closeout(mcdc, idx_work, N_prog, data)

if mcdc["technique"]["domain_decomposition"]:
Expand Down Expand Up @@ -563,7 +552,7 @@ def step_particle(P_arr, data, prog):


@njit
def generate_precursor_particle(DNP_arr, particle_idx, seed_work, prog):
def generate_precursor_particle(DNP_arr, particle_idx, seed, prog):
mcdc = adapt.mcdc_global(prog)
DNP = DNP_arr[0]

Expand All @@ -574,7 +563,7 @@ def generate_precursor_particle(DNP_arr, particle_idx, seed_work, prog):
# Create new particle
P_new_arr = adapt.local_array(1, type_.particle)
P_new = P_new_arr[0]
part_seed = kernel.split_seed(particle_idx, seed_work)
part_seed = kernel.split_seed(particle_idx, seed)
P_new["rng_seed"] = part_seed
P_new["alive"] = True
P_new["w"] = 1.0
Expand All @@ -583,15 +572,21 @@ def generate_precursor_particle(DNP_arr, particle_idx, seed_work, prog):
P_new["x"] = DNP["x"]
P_new["y"] = DNP["y"]
P_new["z"] = DNP["z"]
P_new["t"] = 0.0

# Sample direction
P_new["ux"], P_new["uy"], P_new["uz"] = kernel.sample_isotropic_direction(P_new_arr)

# Get cell and material
P_new["cell_ID"] = -1
P_new["material_ID"] = -1
_ = geometry.locate_particle(P_new_arr, mcdc)

# Get material
_ = geometry.inspect_geometry(P_new_arr, mcdc)
if P_new["cell_ID"] > -1:
material_ID = P_new["material_ID"]
# Skip if particle is lost
else:
if P_new["material_ID"] == -1:
return

material_ID = P_new["material_ID"]
material = mcdc["materials"][material_ID]
G = material["G"]

Expand Down Expand Up @@ -624,9 +619,6 @@ def generate_precursor_particle(DNP_arr, particle_idx, seed_work, prog):

# Accept if it is inside current census index
if P_new["t"] < mcdc["setting"]["census_time"][idx_census]:
# Reduce precursor weight
DNP["w"] -= 1.0

# Skip if it's beyond time boundary
if P_new["t"] > mcdc["setting"]["time_boundary"]:
return
Expand All @@ -640,13 +632,6 @@ def generate_precursor_particle(DNP_arr, particle_idx, seed_work, prog):
break
P_new["g"] = g_out

# Sample direction
(
P_new["ux"],
P_new["uy"],
P_new["uz"],
) = kernel.sample_isotropic_direction(P_new_arr)

# Push to active bank
adapt.add_active(P_new_arr, prog)

Expand All @@ -669,52 +654,33 @@ def source_precursor_closeout(prog, idx_work, N_prog, data):

@njit
def loop_source_precursor(seed, data, mcdc):
# TODO: censussed neutrons seeding is still not reproducible

# Progress bar indicator
N_prog = 0

# =========================================================================
# Sync. RNG skip ahead for reproducibility
# =========================================================================

# Exscan upper estimate of number of particles generated locally
idx_start, N_local, N_global = kernel.bank_scanning_DNP(
mcdc["bank_precursor"], mcdc
)
# TODO: Domain decomposition (see loop_source)

# =========================================================================
# Loop over precursor sources
# =========================================================================
work_start = mcdc["mpi_work_start_precursor"]
work_size = mcdc["mpi_work_size_precursor"]
work_end = work_start + work_size

for idx_work in range(mcdc["mpi_work_size_precursor"]):
for idx_work in range(work_size):
# Get precursor
DNP_arr = mcdc["bank_precursor"]["precursors"][idx_work : (idx_work + 1)]
DNP = DNP_arr[0]
seed_precursor = kernel.split_seed(work_start + idx_work, seed)
# Note the seed is only used once

# Determine number of particles to be generated
w = DNP["w"]
N = math.floor(w)
# "Roulette" the last particle
seed_work = kernel.split_seed(idx_work, seed)
if kernel.rng_from_seed(seed_work) < w - N:
N += 1
DNP["w"] = N
nu = DNP["w"]
xi = kernel.rng_from_seed(seed_precursor)
N = int(math.floor(nu + xi))

# =====================================================================
# Loop over source particles from the source precursor
# =====================================================================

for particle_idx in range(N):

generate_precursor_particle(DNP_arr, particle_idx, seed_work, mcdc)

generate_precursor_particle(DNP_arr, particle_idx, seed_precursor, mcdc)
exhaust_active_bank(data, mcdc)

# =====================================================================
# Closeout
# =====================================================================

source_precursor_closeout(mcdc, idx_work, N_prog, data)


Expand Down
3 changes: 2 additions & 1 deletion mcdc/src/geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,7 +406,8 @@ def report_lost(particle_container):
x = particle["x"]
y = particle["y"]
z = particle["z"]
print("A particle is lost at (", x, y, z, ")")
t = particle["t"]
print("A particle is lost at (", x, y, z, t, ")")
particle["alive"] = False


Expand Down

0 comments on commit 4d5293a

Please sign in to comment.