diff --git a/.gitignore b/.gitignore index a53a3ea6..2da634a5 100644 --- a/.gitignore +++ b/.gitignore @@ -51,4 +51,3 @@ pytestdebug.log docs/build docs/source/pythonapi/generated/ -*.csv diff --git a/examples/fixed_source/slab_reed_dd/input.py b/examples/fixed_source/slab_reed_dd/input.py deleted file mode 100644 index 167e102b..00000000 --- a/examples/fixed_source/slab_reed_dd/input.py +++ /dev/null @@ -1,51 +0,0 @@ -import numpy as np - -import mcdc - -# ============================================================================= -# Set model -# ============================================================================= -# Three slab layers with different materials -# Based on William H. Reed, NSE (1971), 46:2, 309-314, DOI: 10.13182/NSE46-309 - -# Set materials -m1 = mcdc.material(capture=np.array([50.0])) -m2 = mcdc.material(capture=np.array([5.0])) -m3 = mcdc.material(capture=np.array([0.0])) # Vacuum -m4 = mcdc.material(capture=np.array([0.1]), scatter=np.array([[0.9]])) - -# Set surfaces -s1 = mcdc.surface("plane-z", z=0.0, bc="reflective") -s2 = mcdc.surface("plane-z", z=2.0) -s3 = mcdc.surface("plane-z", z=3.0) -s4 = mcdc.surface("plane-z", z=5.0) -s5 = mcdc.surface("plane-z", z=8.0, bc="vacuum") - -# Set cells -mcdc.cell([+s1, -s2], m1) -mcdc.cell([+s2, -s3], m2) -mcdc.cell([+s3, -s4], m3) -mcdc.cell([+s4, -s5], m4) - -# ============================================================================= -# Set source -# ============================================================================= - -# Isotropic source in the absorbing medium -mcdc.source(z=[0.0, 2.0], isotropic=True, prob=50.0) - -# Isotropic source in the first half of the outermost medium, -# with 1/100 strength -mcdc.source(z=[5.0, 6.0], isotropic=True, prob=0.5) - -# ============================================================================= -# Set tally, setting, and run mcdc -# ============================================================================= - -mcdc.tally(scores=["flux"], z=np.linspace(0.0, 8.0, 81)) - -# Setting -mcdc.setting(N_particle=5000) -mcdc.domain_decomposition(z=np.linspace(0.0, 8.0, 5)) -# Run -mcdc.run() diff --git a/examples/fixed_source/slab_reed_dd/process.py b/examples/fixed_source/slab_reed_dd/process.py deleted file mode 100644 index 09799cc3..00000000 --- a/examples/fixed_source/slab_reed_dd/process.py +++ /dev/null @@ -1,123 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import h5py -from scipy.integrate import quad - -# ============================================================================= -# Reference solution (not accurate enough for N_hist > 1E7) -# ============================================================================= - - -def phi1(x): - return ( - 1.0 - - 5.96168047527760 * 10 ** (-47) * np.cosh(52.06761235859028 * x) - - 6.78355315350872 * 10 ** (-56) * np.cosh(62.76152118553390 * x) - - 7.20274049646598 * 10 ** (-84) * np.cosh(95.14161078659372 * x) - - 6.34541150517664 * 10 ** (-238) * np.cosh(272.5766481169758 * x) - ) - - -def phi2(x): - return ( - 1.685808767651539 * 10**3 * np.exp(-5.206761235859028 * x) - + 3.143867366942945 * 10**4 * np.exp(-6.276152118553390 * x) - + 2.879977113018352 * 10**7 * np.exp(-9.514161078659372 * x) - + 8.594190506002560 * 10**22 * np.exp(-27.25766481169758 * x) - + 1.298426035202193 * 10 ** (-36) * np.exp(27.25766481169758 * x) - + 1.432344656303454 * 10 ** (-13) * np.exp(9.514161078659372 * x) - + 1.514562265056083 * 10 ** (-9) * np.exp(6.276152118553390 * x) - + 1.594431209450755 * 10 ** (-8) * np.exp(5.206761235859028 * x) - ) - - -def phi3(x): - return 1.105109108062394 - - -def phi4(x): - return ( - 10.0 - - 0.1983746883968300 * np.exp(0.5254295183311557 * x) - - 7.824765332896027 * 10 ** (-5) * np.exp(1.108937229227813 * x) - - 9.746660212187006 * 10 ** (-6) * np.exp(1.615640334315550 * x) - - 2.895098351422132 * 10 ** (-13) * np.exp(4.554850586269065 * x) - - 75.34793864805979 * np.exp(-0.5254295183311557 * x) - - 20.42874998426011 * np.exp(-1.108937229227813 * x) - - 7.129175418204712 * 10 ** (2) * np.exp(-1.615640334315550 * x) - - 2.716409367577795 * 10 ** (9) * np.exp(-4.554850586269065 * x) - ) - - -def phi5(x): - return ( - 31.53212162577067 * np.exp(-0.5254295183311557 * x) - + 26.25911060454856 * np.exp(-1.108937229227813 * x) - + 1.841223066417334 * 10 ** (3) * np.exp(-1.615640334315550 * x) - + 1.555593549394869 * 10 ** (11) * np.exp(-4.554850586269065 * x) - - 3.119310353653182 * 10 ** (-3) * np.exp(0.5254295183311557 * x) - - 6.336401143340483 * 10 ** (-7) * np.exp(1.108937229227813 * x) - - 3.528757679361232 * 10 ** (-8) * np.exp(1.615640334315550 * x) - - 4.405514335746888 * 10 ** (-18) * np.exp(4.554850586269065 * x) - ) - - -def f_phi(x1, x2): - dx = x2 - x1 - if x2 <= 2.0: - return quad(phi1, x1, x2)[0] / dx - if x2 <= 3.0: - return quad(phi2, x1, x2)[0] / dx - if x2 <= 5.0: - return quad(phi3, x1, x2)[0] / dx - if x2 <= 6.0: - return quad(phi4, x1, x2)[0] / dx - return quad(phi5, x1, x2)[0] / dx - - -def f_phi_x(x): - if x <= 2.0: - return phi1(x) - if x <= 3.0: - return phi2(x) - if x <= 5.0: - return phi3(x) - if x <= 6.0: - return phi4(x) - return phi5(x) - - -with h5py.File("output.h5", "r") as f: - x = f["tally/grid/z"][:] - dx = x[1] - x[0] - x_mid = 0.5 * (x[:-1] + x[1:]) - -phi_ref = np.zeros_like(x_mid) -phi_x_ref = np.zeros_like(x) - -for i in range(len(x)): - phi_x_ref[i] = f_phi_x(x[i]) - -for i in range(len(x_mid)): - phi_ref[i] = f_phi(x[i], x[i + 1]) - -# ============================================================================= -# Plot -# ============================================================================= - -# Load output -with h5py.File("output.h5", "r") as f: - # Note the spatial (dx) and source strength (100+1) normalization - phi = f["tally/flux/mean"][:] / dx * 101.0 - phi_sd = f["tally/flux/sdev"][:] / dx * 101.0 - -# Flux - spatial average -plt.plot(x_mid, phi, "-b", label="MC") -plt.fill_between(x_mid, phi - phi_sd, phi + phi_sd, alpha=0.2, color="b") -plt.plot(x_mid, phi_ref, "--r", label="ref.") -plt.xlabel(r"$x$, cm") -plt.ylabel("Flux") -plt.grid() -plt.legend() -plt.title(r"$\bar{\phi}_i$") -plt.show() diff --git a/mcdc/__init__.py b/mcdc/__init__.py index 632d644a..09595dfb 100644 --- a/mcdc/__init__.py +++ b/mcdc/__init__.py @@ -23,7 +23,6 @@ uq, print_card, reset_cards, - domain_decomposition, ) from mcdc.main import run, prepare from mcdc.visualizer import visualize diff --git a/mcdc/card.py b/mcdc/card.py index bec92980..3e1dcfd1 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -80,12 +80,6 @@ def reset(self): "ww": np.ones([1, 1, 1, 1]), "ww_width": 2.5, "ww_mesh": make_card_mesh(), - "domain_decomposition": False, - "dd_idx": 0, - "dd_mesh": make_card_mesh(), - "dd_exchange_rate": 0, - "dd_repro": False, - "dd_work_ratio": np.array([1]), "weight_roulette": False, "wr_threshold": 0.0, "wr_survive": 1.0, diff --git a/mcdc/constant.py b/mcdc/constant.py index 431d123a..d25df2a4 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -14,7 +14,6 @@ EVENT_TIME_BOUNDARY = 1 << 7 EVENT_LATTICE = 1 << 8 EVENT_SURFACE_MOVE = 1 << 9 -EVENT_DOMAIN = 1 << 10 # Gyration raius type GYRATION_RADIUS_ALL = 0 @@ -39,12 +38,6 @@ PREC = 1.0 + 1e-5 # Precision factor to determine if a distance is smaller BANKMAX = 100 # Default maximum active bank -# Domain Decomp mesh crossing flags -MESH_X = 0 -MESH_Y = 1 -MESH_Z = 2 -MESH_T = 3 - # RNG LCG parameters RNG_G = nb.uint64(2806196910506780709) RNG_C = nb.uint64(1) diff --git a/mcdc/input_.py b/mcdc/input_.py index 2bc1fd5d..0ec2022f 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -1173,50 +1173,6 @@ def weight_window(x=None, y=None, z=None, t=None, window=None, width=None): return card -def domain_decomposition( - x=None, - y=None, - z=None, - t=None, - exchange_rate=100, - bank_size=1e5, - work_ratio=None, - repro=True, -): - card = mcdc.input_deck.technique - card["domain_decomposition"] = True - card["domain_bank_size"] = int(1e5) - card["dd_exchange_rate"] = int(exchange_rate) - card["dd_repro"] = repro - dom_num = 1 - # Set mesh - if x is not None: - card["dd_mesh"]["x"] = x - dom_num *= len(x) - if y is not None: - card["dd_mesh"]["y"] = y - dom_num *= len(y) - if z is not None: - card["dd_mesh"]["z"] = z - dom_num += len(z) - if t is not None: - card["dd_mesh"]["t"] = t - dom_num += len(t) - # Set work ratio - if work_ratio is None: - card["dd_work_ratio"] = None - elif work_ratio is not None: - card["dd_work_ratio"] = work_ratio - card["dd_idx"] = 0 - card["dd_xp_neigh"] = [] - card["dd_xn_neigh"] = [] - card["dd_yp_neigh"] = [] - card["dd_yn_neigh"] = [] - card["dd_zp_neigh"] = [] - card["dd_zn_neigh"] = [] - return card - - def iQMC( g=None, t=None, diff --git a/mcdc/kernel.py b/mcdc/kernel.py index f5fb68bb..074dfc5e 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -7,562 +7,11 @@ import mcdc.type_ as type_ from mcdc.constant import * -from mcdc.print_ import print_error, print_msg +from mcdc.print_ import print_error from mcdc.type_ import score_list, iqmc_score_list from mcdc.loop import loop_source -# ============================================================================= -# Domain Decomposition -# ============================================================================= - -# ============================================================================= -# Domain crossing event -# ============================================================================= - - -@njit -def domain_crossing(P, mcdc): - # Domain mesh crossing - max_size = mcdc["technique"]["dd_exchange_rate"] - if mcdc["technique"]["domain_decomposition"]: - mesh = mcdc["technique"]["dd_mesh"] - # Determine which dimension is crossed - x, y, z, t, directions = mesh_crossing_evaluate(P, mesh) - if len(directions) == 0: - return - elif len(directions) > 1: - for direction in directions[1:]: - if direction == MESH_X: - P["x"] -= SHIFT * P["ux"] / np.abs(P["ux"]) - if direction == MESH_Y: - P["y"] -= SHIFT * P["uy"] / np.abs(P["uy"]) - if direction == MESH_Z: - P["z"] -= SHIFT * P["uz"] / np.abs(P["uz"]) - flag = directions[0] - # Score on tally - if flag == MESH_X and P["ux"] > 0: - add_particle(copy_particle(P), mcdc["bank_domain_xp"]) - if mcdc["bank_domain_xp"]["size"] == max_size: - dd_particle_send(mcdc) - if flag == MESH_X and P["ux"] < 0: - add_particle(copy_particle(P), mcdc["bank_domain_xn"]) - if mcdc["bank_domain_xn"]["size"] == max_size: - dd_particle_send(mcdc) - if flag == MESH_Y and P["uy"] > 0: - add_particle(copy_particle(P), mcdc["bank_domain_yp"]) - if mcdc["bank_domain_yp"]["size"] == max_size: - dd_particle_send(mcdc) - if flag == MESH_Y and P["uy"] < 0: - add_particle(copy_particle(P), mcdc["bank_domain_yn"]) - if mcdc["bank_domain_yn"]["size"] == max_size: - dd_particle_send(mcdc) - if flag == MESH_Z and P["uz"] > 0: - add_particle(copy_particle(P), mcdc["bank_domain_zp"]) - if mcdc["bank_domain_zp"]["size"] == max_size: - dd_particle_send(mcdc) - if flag == MESH_Z and P["uz"] < 0: - add_particle(copy_particle(P), mcdc["bank_domain_zn"]) - if mcdc["bank_domain_zn"]["size"] == max_size: - dd_particle_send(mcdc) - P["alive"] = False - - -# ============================================================================= -# Send full domain bank -# ============================================================================= - - -@njit -def dd_particle_send(mcdc): - with objmode(): - for i in range( - max( - len(mcdc["technique"]["dd_xp_neigh"]), - len(mcdc["technique"]["dd_xn_neigh"]), - len(mcdc["technique"]["dd_yp_neigh"]), - len(mcdc["technique"]["dd_yn_neigh"]), - len(mcdc["technique"]["dd_zp_neigh"]), - len(mcdc["technique"]["dd_zn_neigh"]), - ) - ): - if mcdc["technique"]["dd_xp_neigh"].size > i: - size = mcdc["bank_domain_xp"]["size"] - ratio = int(size / len(mcdc["technique"]["dd_xp_neigh"])) - start = ratio * i - end = start + ratio - if i == len(mcdc["technique"]["dd_xp_neigh"]) - 1: - end = size - bank = np.array(mcdc["bank_domain_xp"]["particles"][start:end]) - request1 = MPI.COMM_WORLD.send( - bank, dest=mcdc["technique"]["dd_xp_neigh"][i], tag=1 - ) - - if mcdc["technique"]["dd_xn_neigh"].size > i: - size = mcdc["bank_domain_xn"]["size"] - ratio = int(size / len(mcdc["technique"]["dd_xn_neigh"])) - start = ratio * i - end = start + ratio - if i == len(mcdc["technique"]["dd_xn_neigh"]) - 1: - end = size - bank = np.array(mcdc["bank_domain_xn"]["particles"][start:end]) - request2 = MPI.COMM_WORLD.send( - bank, dest=mcdc["technique"]["dd_xn_neigh"][i], tag=2 - ) - - if mcdc["technique"]["dd_yp_neigh"].size > i: - size = mcdc["bank_domain_yp"]["size"] - ratio = int(size / len(mcdc["technique"]["dd_yp_neigh"])) - start = ratio * i - end = start + ratio - if i == len(mcdc["technique"]["dd_yp_neigh"]) - 1: - end = size - bank = np.array(mcdc["bank_domain_yp"]["particles"][start:end]) - request3 = MPI.COMM_WORLD.send( - bank, dest=mcdc["technique"]["dd_yp_neigh"][i], tag=3 - ) - - if mcdc["technique"]["dd_yn_neigh"].size > i: - size = mcdc["bank_domain_yn"]["size"] - ratio = int(size / len(mcdc["technique"]["dd_yn_neigh"])) - start = ratio * i - end = start + ratio - if i == len(mcdc["technique"]["dd_yn_neigh"]) - 1: - end = size - bank = np.array(mcdc["bank_domain_yn"]["particles"][start:end]) - request4 = MPI.COMM_WORLD.send( - bank, dest=mcdc["technique"]["dd_yn_neigh"][i], tag=4 - ) - - if mcdc["technique"]["dd_zp_neigh"].size > i: - size = mcdc["bank_domain_zp"]["size"] - ratio = int(size / len(mcdc["technique"]["dd_zp_neigh"])) - start = ratio * i - end = start + ratio - if i == len(mcdc["technique"]["dd_zp_neigh"]) - 1: - end = size - bank = np.array(mcdc["bank_domain_zp"]["particles"][start:end]) - request5 = MPI.COMM_WORLD.send( - bank, dest=mcdc["technique"]["dd_zp_neigh"][i], tag=5 - ) - - if mcdc["technique"]["dd_zn_neigh"].size > i: - size = mcdc["bank_domain_zn"]["size"] - ratio = int(size / len(mcdc["technique"]["dd_zn_neigh"])) - start = ratio * i - end = start + ratio - if i == len(mcdc["technique"]["dd_zn_neigh"]) - 1: - end = size - bank = np.array(mcdc["bank_domain_zn"]["particles"][start:end]) - request6 = MPI.COMM_WORLD.send( - bank, dest=mcdc["technique"]["dd_zn_neigh"][i], tag=6 - ) - - sent_particles = ( - mcdc["bank_domain_xp"]["size"] - + mcdc["bank_domain_xn"]["size"] - + mcdc["bank_domain_yp"]["size"] - + mcdc["bank_domain_yn"]["size"] - + mcdc["bank_domain_zp"]["size"] - + mcdc["bank_domain_zn"]["size"] - ) - mcdc["technique"]["dd_sent"] += sent_particles - - mcdc["bank_domain_xp"]["size"] = 0 - mcdc["bank_domain_xn"]["size"] = 0 - mcdc["bank_domain_yp"]["size"] = 0 - mcdc["bank_domain_yn"]["size"] = 0 - mcdc["bank_domain_zp"]["size"] = 0 - mcdc["bank_domain_zn"]["size"] = 0 - - -# ============================================================================= -# Recieve particles and clear banks -# ============================================================================= - - -@njit -def dd_particle_receive(mcdc): - buff = np.zeros( - mcdc["bank_domain_xp"]["particles"].shape[0], dtype=type_.particle_record - ) - - with objmode(size="int64"): - bankr = mcdc["bank_active"]["particles"][:0] - size_old = bankr.shape[0] - for i in range( - max( - len(mcdc["technique"]["dd_xp_neigh"]), - len(mcdc["technique"]["dd_xn_neigh"]), - len(mcdc["technique"]["dd_yp_neigh"]), - len(mcdc["technique"]["dd_yn_neigh"]), - len(mcdc["technique"]["dd_zp_neigh"]), - len(mcdc["technique"]["dd_zn_neigh"]), - ) - ): - if mcdc["technique"]["dd_xp_neigh"].size > i: - received1 = MPI.COMM_WORLD.irecv( - source=mcdc["technique"]["dd_xp_neigh"][i], tag=2 - ) - if received1.Get_status(): - bankr = np.append(bankr, received1.wait()) - else: - MPI.Request.cancel(received1) - - if mcdc["technique"]["dd_xn_neigh"].size > i: - received2 = MPI.COMM_WORLD.irecv( - source=mcdc["technique"]["dd_xn_neigh"][i], tag=1 - ) - if received2.Get_status(): - bankr = np.append(bankr, received2.wait()) - else: - MPI.Request.cancel(received2) - - if mcdc["technique"]["dd_yp_neigh"].size > i: - received3 = MPI.COMM_WORLD.irecv( - source=mcdc["technique"]["dd_yp_neigh"][i], tag=4 - ) - if received3.Get_status(): - bankr = np.append(bankr, received3.wait()) - else: - MPI.Request.cancel(received3) - - if mcdc["technique"]["dd_yn_neigh"].size > i: - received4 = MPI.COMM_WORLD.irecv( - source=mcdc["technique"]["dd_yn_neigh"][i], tag=3 - ) - if received4.Get_status(): - bankr = np.append(bankr, received4.wait()) - else: - MPI.Request.cancel(received4) - - if mcdc["technique"]["dd_zp_neigh"].size > i: - received5 = MPI.COMM_WORLD.irecv( - source=mcdc["technique"]["dd_zp_neigh"][i], tag=6 - ) - if received5.Get_status(): - bankr = np.append(bankr, received5.wait()) - else: - MPI.Request.cancel(received5) - - if mcdc["technique"]["dd_zn_neigh"].size > i: - received6 = MPI.COMM_WORLD.irecv( - source=mcdc["technique"]["dd_zn_neigh"][i], tag=5 - ) - if received6.Get_status(): - bankr = np.append(bankr, received6.wait()) - else: - MPI.Request.cancel(received6) - - size = bankr.shape[0] - # Set output buffer - for i in range(size): - buff[i] = bankr[i] - - # Set source bank from buffer - for i in range(size): - add_particle(buff[i], mcdc["bank_active"]) - mcdc["technique"]["dd_sent"] -= size - - -# ============================================================================= -# Particle in domain -# ============================================================================= - - -# Check if particle is in domain -@njit -def particle_in_domain(P, mcdc): - 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 - - d_iz = int(d_idx / (d_Nx * d_Ny)) - d_iy = int((d_idx - d_Nx * d_Ny * d_iz) / d_Nx) - d_ix = int(d_idx - d_Nx * d_Ny * d_iz - d_Nx * d_iy) - - x_cell = binary_search(P["x"], mcdc["technique"]["dd_mesh"]["x"]) - y_cell = binary_search(P["y"], mcdc["technique"]["dd_mesh"]["y"]) - z_cell = binary_search(P["z"], mcdc["technique"]["dd_mesh"]["z"]) - - if d_ix == x_cell: - if d_iy == y_cell: - if d_iz == z_cell: - return True - return False - - -# ============================================================================= -# Source in domain -# ============================================================================= - - -# Check for source in domain -@njit -def source_in_domain(source, domain_mesh, d_idx): - d_Nx = domain_mesh["x"].size - 1 - d_Ny = domain_mesh["y"].size - 1 - d_Nz = domain_mesh["z"].size - 1 - - d_iz = int(d_idx / (d_Nx * d_Ny)) - d_iy = int((d_idx - d_Nx * d_Ny * d_iz) / d_Nx) - d_ix = int(d_idx - d_Nx * d_Ny * d_iz - d_Nx * d_iy) - - d_x = [domain_mesh["x"][d_ix], domain_mesh["x"][d_ix + 1]] - d_y = [domain_mesh["y"][d_iy], domain_mesh["y"][d_iy + 1]] - d_z = [domain_mesh["z"][d_iz], domain_mesh["z"][d_iz + 1]] - - if ( - d_x[0] <= source["box_x"][0] <= d_x[1] - or d_x[0] <= source["box_x"][1] <= d_x[1] - or (source["box_x"][0] < d_x[0] and source["box_x"][1] > d_x[1]) - ): - if ( - d_y[0] <= source["box_y"][0] <= d_y[1] - or d_y[0] <= source["box_y"][1] <= d_y[1] - or (source["box_y"][0] < d_y[0] and source["box_y"][1] > d_y[1]) - ): - if ( - d_z[0] <= source["box_z"][0] <= d_z[1] - or d_z[0] <= source["box_z"][1] <= d_z[1] - or (source["box_z"][0] < d_z[0] and source["box_z"][1] > d_z[1]) - ): - return True - else: - return False - else: - return False - else: - return False - - -# ============================================================================= -# Compute domain load -# ============================================================================= - - -@njit -def domain_work(mcdc, domain, N): - domain_mesh = mcdc["technique"]["dd_mesh"] - - d_Nx = domain_mesh["x"].size - 1 - d_Ny = domain_mesh["y"].size - 1 - d_Nz = domain_mesh["z"].size - 1 - work_start = 0 - for d_idx in range(domain): - d_iz = int(d_idx / (d_Nx * d_Ny)) - d_iy = int((d_idx - d_Nx * d_Ny * d_iz) / d_Nx) - d_ix = int(d_idx - d_Nx * d_Ny * d_iz - d_Nx * d_iy) - - d_x = [domain_mesh["x"][d_ix], domain_mesh["x"][d_ix + 1]] - d_y = [domain_mesh["y"][d_iy], domain_mesh["y"][d_iy + 1]] - d_z = [domain_mesh["z"][d_iz], domain_mesh["z"][d_iz + 1]] - # Compute volumes of sources and numbers of particles - - Psum = 0 - - Nm = 0 - num_source = 0 - for source in mcdc["sources"]: - Psum += source["prob"] - num_source += 1 - Vi = np.zeros(num_source) - Vim = np.zeros(num_source) - Ni = np.zeros(num_source) - i = 0 - for source in mcdc["sources"]: - Ni[i] = N * source["prob"] / Psum - Vi[i] = 1 - Vim[i] = 1 - if source["box"] == True: - xV = source["box_x"][1] - source["box_x"][0] - if xV != 0: - Vi[i] *= xV - Vim[i] *= min(source["box_x"][1], d_x[1]) - max( - source["box_x"][0], d_x[0] - ) - yV = source["box_y"][1] - source["box_y"][0] - if yV != 0: - Vi[i] *= yV - Vim[i] *= min(source["box_y"][1], d_y[1]) - max( - source["box_y"][0], d_y[0] - ) - zV = source["box_z"][1] - source["box_z"][0] - if zV != 0: - Vi[i] *= zV - Vim[i] *= min(source["box_z"][1], d_z[1]) - max( - source["box_z"][0], d_z[0] - ) - if not source_in_domain(source, domain_mesh, d_idx): - Vim[i] = 0 - i += 1 - for source in range(num_source): - Nm += Ni[source] * Vim[source] / Vi[source] - work_start += Nm - d_idx = domain - d_iz = int(mcdc["dd_idx"] / (d_Nx * d_Ny)) - d_iy = int((mcdc["dd_idx"] - d_Nx * d_Ny * d_iz) / d_Nx) - d_ix = int(mcdc["dd_idx"] - d_Nx * d_Ny * d_iz - d_Nx * d_iy) - - d_x = [domain_mesh["x"][d_ix], domain_mesh["x"][d_ix + 1]] - d_y = [domain_mesh["y"][d_iy], domain_mesh["y"][d_iy + 1]] - d_z = [domain_mesh["z"][d_iz], domain_mesh["z"][d_iz + 1]] - # Compute volumes of sources and numbers of particles - num_source = len(mcdc["sources"]) - Vi = np.zeros(num_source) - Vim = np.zeros(num_source) - Ni = np.zeros(num_source) - Psum = 0 - - Nm = 0 - for source in mcdc["sources"]: - Psum += source["prob"] - i = 0 - for source in mcdc["sources"]: - Ni[i] = N * source["prob"] / Psum - Vi[i] = 1 - Vim[i] = 1 - - if source["box"] == True: - xV = source["box_x"][1] - source["box_x"][0] - if xV != 0: - Vi[i] *= xV - Vim[i] *= min(source["box_x"][1], d_x[1]) - max( - source["box_x"][0], d_x[0] - ) - yV = source["box_y"][1] - source["box_y"][0] - if yV != 0: - Vi[i] *= yV - Vim[i] *= min(source["box_y"][1], d_y[1]) - max( - source["box_y"][0], d_y[0] - ) - zV = source["box_z"][1] - source["box_z"][0] - if zV != 0: - Vi[i] *= zV - Vim[i] *= min(source["box_z"][1], d_z[1]) - max( - source["box_z"][0], d_z[0] - ) - i += 1 - for source in range(num_source): - Nm += Ni[source] * Vim[source] / Vi[source] - Nm /= mcdc["technique"]["dd_work_ratio"][domain] - rank = mcdc["mpi_rank"] - if mcdc["technique"]["dd_work_ratio"][domain] > 1: - work_start += Nm * (rank - np.sum(mcdc["technique"]["dd_work_ratio"][0:d_idx])) - total_v = 0 - for source in range(len(mcdc["sources"])): - total_v += Vim[source] - i = 0 - for source in mcdc["sources"]: - if total_v != 0: - source["prob"] *= 2 * Vim[i] / total_v - i += 1 - return (int(Nm), int(work_start)) - - -# ============================================================================= -# Source particle in domain only -# ============================================================================= - - -@njit() -def source_particle_dd(seed, mcdc): - domain_mesh = mcdc["technique"]["dd_mesh"] - d_idx = mcdc["dd_idx"] - - d_Nx = domain_mesh["x"].size - 1 - d_Ny = domain_mesh["y"].size - 1 - d_Nz = domain_mesh["z"].size - 1 - - d_iz = int(mcdc["dd_idx"] / (d_Nx * d_Ny)) - d_iy = int((mcdc["dd_idx"] - d_Nx * d_Ny * d_iz) / d_Nx) - d_ix = int(mcdc["dd_idx"] - d_Nx * d_Ny * d_iz - d_Nx * d_iy) - - d_x = [domain_mesh["x"][d_ix], domain_mesh["x"][d_ix + 1]] - d_y = [domain_mesh["y"][d_iy], domain_mesh["y"][d_iy + 1]] - d_z = [domain_mesh["z"][d_iz], domain_mesh["z"][d_iz + 1]] - - P = np.zeros(1, dtype=type_.particle_record)[0] - - P["rng_seed"] = seed - # Sample source - xi = rng(P) - tot = 0.0 - for source in mcdc["sources"]: - if source_in_domain(source, domain_mesh, d_idx): - tot += source["prob"] - if tot >= xi: - break - - # Position - if source["box"]: - x = sample_uniform( - max(source["box_x"][0], d_x[0]), min(source["box_x"][1], d_x[1]), P - ) - y = sample_uniform( - max(source["box_y"][0], d_y[0]), min(source["box_y"][1], d_y[1]), P - ) - z = sample_uniform( - max(source["box_z"][0], d_z[0]), min(source["box_z"][1], d_z[1]), P - ) - - else: - x = source["x"] - y = source["y"] - z = source["z"] - - # Direction - if source["isotropic"]: - ux, uy, uz = sample_isotropic_direction(P) - elif source["white"]: - ux, uy, uz = sample_white_direction( - source["white_x"], source["white_y"], source["white_z"], P - ) - else: - ux = source["ux"] - uy = source["uy"] - uz = source["uz"] - - # Energy and time - g = sample_discrete(source["group"], P) - t = sample_uniform(source["time"][0], source["time"][1], P) - - # Make and return particle - P["x"] = x - P["y"] = y - P["z"] = z - P["t"] = t - P["ux"] = ux - P["uy"] = uy - P["uz"] = uz - P["g"] = g - P["w"] = 1 - P["sensitivity_ID"] = 0 - return P - - -@njit -def distribute_work_dd(N, mcdc, precursor=False): - # Total # of work - work_size_total = N - - if not mcdc["technique"]["dd_repro"]: - work_size, work_start = domain_work(mcdc, mcdc["dd_idx"], N) - else: - work_start = 0 - work_size = work_size_total - - if not precursor: - mcdc["mpi_work_start"] = work_start - mcdc["mpi_work_size"] = work_size - mcdc["mpi_work_size_total"] = work_size_total - else: - mcdc["mpi_work_start_precursor"] = work_start - mcdc["mpi_work_size_precursor"] = work_size - mcdc["mpi_work_size_total_precursor"] = work_size_total - - # ============================================================================= # Random sampling # ============================================================================= @@ -880,8 +329,7 @@ def manage_particle_banks(seed, mcdc): ] # MPI rebalance - if not mcdc["technique"]["domain_decomposition"]: - bank_rebalance(mcdc) + bank_rebalance(mcdc) # Zero out census bank mcdc["bank_census"]["size"] = 0 @@ -1044,14 +492,6 @@ def total_weight(bank): return buff[0] -@njit -def allreduce(value): - total = np.zeros(1, np.float64) - with objmode(): - MPI.COMM_WORLD.Allreduce(np.array([value], np.float64), total, MPI.SUM) - return total[0] - - @njit def bank_rebalance(mcdc): # Scan the bank @@ -1845,32 +1285,6 @@ def mesh_uniform_get_index(P, mesh, trans): return x, y, z -@njit -def mesh_crossing_evaluate(P, mesh): - # Shift backward - shift_particle(P, -2 * SHIFT) - t1, x1, y1, z1, outside1 = mesh_get_index(P, mesh) - - # Double shift forward - shift_particle(P, 4 * SHIFT) - t2, x2, y2, z2, outside2 = mesh_get_index(P, mesh) - - # Return particle to initial position - shift_particle(P, -2 * SHIFT) - - # Determine dimension crossed - directions = [] - - if x1 != x2: - directions.append(MESH_X) - if y1 != y2: - directions.append(MESH_Y) - if z1 != z2: - directions.append(MESH_Z) - - return x1, y1, z1, t1, directions - - # ============================================================================= # Tally operations # ============================================================================= @@ -2261,10 +1675,6 @@ def move_to_event(P, mcdc): if mcdc["cycle_active"]: d_mesh = distance_to_mesh(P, mcdc["tally"]["mesh"], mcdc) - d_domain = INF - if mcdc["cycle_active"] and mcdc["technique"]["domain_decomposition"]: - d_domain = distance_to_mesh(P, mcdc["technique"]["dd_mesh"], mcdc) - if mcdc["technique"]["iQMC"]: d_iqmc_mesh = distance_to_mesh(P, mcdc["technique"]["iqmc"]["mesh"], mcdc) if d_iqmc_mesh < d_mesh: @@ -2291,9 +1701,7 @@ def move_to_event(P, mcdc): # ========================================================================= # Find the minimum - distance = min( - d_boundary, d_time_boundary, d_time_census, d_mesh, d_collision, d_domain - ) + distance = min(d_boundary, d_time_boundary, d_time_census, d_mesh, d_collision) # Remove the boundary event if it is not the nearest if d_boundary > distance * PREC: event = 0 @@ -2305,8 +1713,6 @@ def move_to_event(P, mcdc): event += EVENT_CENSUS if d_mesh <= distance * PREC: event += EVENT_MESH - if d_domain <= distance * PREC: - event += EVENT_DOMAIN if d_collision == distance: event = EVENT_COLLISION diff --git a/mcdc/loop.py b/mcdc/loop.py index e550a38c..202edf34 100644 --- a/mcdc/loop.py +++ b/mcdc/loop.py @@ -165,22 +165,10 @@ def loop_source(seed, mcdc): # Check if it is beyond current census index idx_census = mcdc["idx_census"] if P["t"] > mcdc["setting"]["census_time"][idx_census]: - if mcdc["technique"]["domain_decomposition"]: - if mcdc["technique"]["dd_work_ratio"][mcdc["dd_idx"]] > 0: - P["w"] /= mcdc["technique"]["dd_work_ratio"][mcdc["dd_idx"]] - if kernel.particle_in_domain(P, mcdc): - kernel.add_particle(P, mcdc["bank_census"]) - else: - kernel.add_particle(P, mcdc["bank_census"]) + kernel.add_particle(P, mcdc["bank_census"]) else: # Add the source particle into the active bank - if mcdc["technique"]["domain_decomposition"]: - if mcdc["technique"]["dd_work_ratio"][mcdc["dd_idx"]] > 0: - P["w"] /= mcdc["technique"]["dd_work_ratio"][mcdc["dd_idx"]] - if kernel.particle_in_domain(P, mcdc): - kernel.add_particle(P, mcdc["bank_active"]) - else: - kernel.add_particle(P, mcdc["bank_active"]) + kernel.add_particle(P, mcdc["bank_active"]) # ===================================================================== # Run the source particle and its secondaries @@ -221,59 +209,6 @@ def loop_source(seed, mcdc): with objmode(): print_progress(percent, mcdc) - if mcdc["technique"]["domain_decomposition"]: - kernel.dd_particle_send(mcdc) - terminated = False - max_work = 1 - kernel.dd_particle_receive(mcdc) - while not terminated: - if mcdc["bank_active"]["size"] > 0: - # Loop until active bank is exhausted - while mcdc["bank_active"]["size"] > 0: - - P = kernel.get_particle(mcdc["bank_active"], mcdc) - if not kernel.particle_in_domain(P, mcdc) and P["alive"] == True: - print("recieved particle not in domain, position:") - - # Apply weight window - if mcdc["technique"]["weight_window"]: - kernel.weight_window(P, mcdc) - - # Particle tracker - if mcdc["setting"]["track_particle"]: - mcdc["particle_track_particle_ID"] += 1 - - # Particle loop - loop_particle(P, mcdc) - - # Tally history closeout for one-batch fixed-source simulation - if ( - not mcdc["setting"]["mode_eigenvalue"] - and mcdc["setting"]["N_batch"] == 1 - ): - kernel.tally_closeout_history(mcdc) - - # Send all domain particle banks - kernel.dd_particle_send(mcdc) - - # Check for incoming particles - kernel.dd_particle_receive(mcdc) - work_remaining = int(kernel.allreduce(mcdc["bank_active"]["size"])) - total_sent = int(kernel.allreduce(mcdc["technique"]["dd_sent"])) - if work_remaining > max_work: - max_work = work_remaining - - # Progress printout - """ - percent = 1 - work_remaining / max_work - if mcdc["setting"]["progress_bar"] and int(percent * 100.0) > N_prog: - N_prog += 1 - with objmode(): - print_progress(percent, mcdc) - """ - if work_remaining + total_sent == 0: - terminated = True - # ========================================================================= # Particle loop @@ -333,18 +268,10 @@ def loop_particle(P, mcdc): # Surface crossing if event & EVENT_SURFACE: kernel.surface_crossing(P, mcdc) - if event & EVENT_DOMAIN: - if not ( - mcdc["surfaces"][P["surface_ID"]]["reflective"] - or mcdc["surfaces"][P["surface_ID"]]["vacuum"] - ): - kernel.domain_crossing(P, mcdc) # Lattice or mesh crossing (skipped if surface crossing) elif event & EVENT_LATTICE or event & EVENT_MESH: kernel.shift_particle(P, SHIFT) - if event & EVENT_DOMAIN: - kernel.domain_crossing(P, mcdc) # Moving surface transition if event & EVENT_SURFACE_MOVE: diff --git a/mcdc/main.py b/mcdc/main.py index 1533efc1..cfe5eead 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -137,127 +137,6 @@ def run(): closeout(mcdc) -# ============================================================================= -# prepare domain decomposition -# ============================================================================= -def get_d_idx(i, j, k, ni, nj): - N = i + j * ni + k * ni * nj - return N - - -def get_indexes(N, nx, ny): - k = int(N / (nx * ny)) - j = int((N - nx * ny * k) / nx) - i = int(N - nx * ny * k - nx * j) - return i, j, k - - -def get_neighbors(N, w, nx, ny, nz): - i, j, k = get_indexes(N, nx, ny) - if i > 0: - xn = get_d_idx(i - 1, j, k, nx, ny) - else: - xn = None - if i < (nx - 1): - xp = get_d_idx(i + 1, j, k, nx, ny) - else: - xp = None - if j > 0: - yn = get_d_idx(i, j - 1, k, nx, ny) - else: - yn = None - if j < (ny - 1): - yp = get_d_idx(i, j + 1, k, nx, ny) - else: - yp = None - if k > 0: - zn = get_d_idx(i, j, k - 1, nx, ny) - else: - zn = None - if k < (nz - 1): - zp = get_d_idx(i, j, k + 1, nx, ny) - else: - zp = None - return xn, xp, yn, yp, zn, zp - - -def dd_prepare(): - work_ratio = input_deck.technique["dd_work_ratio"] - - d_Nx = input_deck.technique["dd_mesh"]["x"].size - 1 - d_Ny = input_deck.technique["dd_mesh"]["y"].size - 1 - d_Nz = input_deck.technique["dd_mesh"]["z"].size - 1 - - input_deck.setting["bank_active_buff"] = 1000 - if input_deck.technique["dd_exchange_rate"] == None: - input_deck.technique["dd_exchange_rate"] = 100 - - if work_ratio is None: - work_ratio = np.ones(d_Nx * d_Ny * d_Nz) - input_deck.technique["dd_work_ratio"] = work_ratio - - if ( - input_deck.technique["domain_decomposition"] - and np.sum(work_ratio) != MPI.COMM_WORLD.Get_size() - ): - print_msg( - "Domain work ratio not equal to number of processors, %i != %i " - % (np.sum(work_ratio), MPI.COMM_WORLD.Get_size()) - ) - exit() - - if input_deck.technique["domain_decomposition"]: - # Assigning domain index - i = 0 - rank_info = [] - for n in range(d_Nx * d_Ny * d_Nz): - ranks = [] - for r in range(int(work_ratio[n])): - ranks.append(i) - if MPI.COMM_WORLD.Get_rank() == i: - d_idx = n - i += 1 - rank_info.append(ranks) - input_deck.technique["dd_idx"] = d_idx - xn, xp, yn, yp, zn, zp = get_neighbors(d_idx, 0, d_Nx, d_Ny, d_Nz) - else: - input_deck.technique["dd_idx"] = 0 - input_deck.technique["dd_xp_neigh"] = [] - input_deck.technique["dd_xn_neigh"] = [] - input_deck.technique["dd_yp_neigh"] = [] - input_deck.technique["dd_yn_neigh"] = [] - input_deck.technique["dd_zp_neigh"] = [] - input_deck.technique["dd_zn_neigh"] = [] - return - - if xp is not None: - input_deck.technique["dd_xp_neigh"] = rank_info[xp] - else: - input_deck.technique["dd_xp_neigh"] = [] - if xn is not None: - input_deck.technique["dd_xn_neigh"] = rank_info[xn] - else: - input_deck.technique["dd_xn_neigh"] = [] - - if yp is not None: - input_deck.technique["dd_yp_neigh"] = rank_info[yp] - else: - input_deck.technique["dd_yp_neigh"] = [] - if yn is not None: - input_deck.technique["dd_yn_neigh"] = rank_info[yn] - else: - input_deck.technique["dd_yn_neigh"] = [] - - if zp is not None: - input_deck.technique["dd_zp_neigh"] = rank_info[zp] - else: - input_deck.technique["dd_zp_neigh"] = [] - if zn is not None: - input_deck.technique["dd_zn_neigh"] = rank_info[zn] - else: - input_deck.technique["dd_zn_neigh"] = [] - - def prepare(): """ Preparing the MC transport simulation: @@ -266,8 +145,6 @@ def prepare(): (3) Create and set up global variable container `mcdc` """ - dd_prepare() - # ========================================================================= # Adapt kernels # ========================================================================= @@ -519,7 +396,6 @@ def prepare(): "implicit_capture", "population_control", "weight_window", - "domain_decomposition", "weight_roulette", "iQMC", "IC_generator", @@ -571,31 +447,6 @@ def prepare(): # Survival probability mcdc["technique"]["wr_survive"] = input_deck.technique["wr_survive"] - # ========================================================================= - # Domain Decomposition - # ========================================================================= - - # Set domain mesh - if input_deck.technique["domain_decomposition"]: - name = "dd_mesh" - mcdc["technique"][name]["x"] = input_deck.technique[name]["x"] - mcdc["technique"][name]["y"] = input_deck.technique[name]["y"] - mcdc["technique"][name]["z"] = input_deck.technique[name]["z"] - mcdc["technique"][name]["t"] = input_deck.technique[name]["t"] - mcdc["technique"][name]["mu"] = input_deck.technique[name]["mu"] - mcdc["technique"][name]["azi"] = input_deck.technique[name]["azi"] - # Set exchange rate - mcdc["technique"]["dd_exchange_rate"] = input_deck.technique["dd_exchange_rate"] - mcdc["technique"]["dd_repro"] = input_deck.technique["dd_repro"] - # Set domain index - mcdc["dd_idx"] = input_deck.technique["dd_idx"] - mcdc["technique"]["dd_xp_neigh"] = input_deck.technique["dd_xp_neigh"] - mcdc["technique"]["dd_xn_neigh"] = input_deck.technique["dd_xn_neigh"] - mcdc["technique"]["dd_yp_neigh"] = input_deck.technique["dd_yp_neigh"] - mcdc["technique"]["dd_yn_neigh"] = input_deck.technique["dd_yn_neigh"] - mcdc["technique"]["dd_zp_neigh"] = input_deck.technique["dd_zp_neigh"] - mcdc["technique"]["dd_zn_neigh"] = input_deck.technique["dd_zn_neigh"] - mcdc["technique"]["dd_work_ratio"] = input_deck.technique["dd_work_ratio"] # ========================================================================= # Quasi Monte Carlo @@ -736,10 +587,7 @@ def prepare(): mcdc["mpi_master"] = mcdc["mpi_rank"] == 0 # Distribute work to MPI ranks - if mcdc["technique"]["domain_decomposition"]: - kernel.distribute_work_dd(mcdc["setting"]["N_particle"], mcdc) - else: - kernel.distribute_work(mcdc["setting"]["N_particle"], mcdc) + kernel.distribute_work(mcdc["setting"]["N_particle"], mcdc) # ========================================================================= # Particle banks diff --git a/mcdc/type_.py b/mcdc/type_.py index 7aa92bcf..db10b516 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -694,7 +694,6 @@ def make_type_technique(input_deck): ("iQMC", bool_), ("IC_generator", bool_), ("branchless_collision", bool_), - ("domain_decomposition", bool_), ("uq", bool_), ] @@ -704,24 +703,6 @@ def make_type_technique(input_deck): struct += [("pct", int64), ("pc_factor", float64)] - # ========================================================================= - # domain decomp - # ========================================================================= - # Mesh - mesh, Nx, Ny, Nz, Nt, Nmu, N_azi, Ng = make_type_mesh(card["dd_mesh"]) - struct += [("dd_mesh", mesh)] - struct += [("dd_idx", int64)] - struct += [("dd_sent", int64)] - struct += [("dd_work_ratio", int64, (len(card["dd_work_ratio"]),))] - struct += [("dd_exchange_rate", int64)] - struct += [("dd_repro", bool_)] - struct += [("dd_xp_neigh", int64, (len(card["dd_xp_neigh"]),))] - struct += [("dd_xn_neigh", int64, (len(card["dd_xn_neigh"]),))] - struct += [("dd_yp_neigh", int64, (len(card["dd_yp_neigh"]),))] - struct += [("dd_yn_neigh", int64, (len(card["dd_yn_neigh"]),))] - struct += [("dd_zp_neigh", int64, (len(card["dd_zp_neigh"]),))] - struct += [("dd_zn_neigh", int64, (len(card["dd_zn_neigh"]),))] - # ========================================================================= # Weight window # ========================================================================= @@ -1075,22 +1056,6 @@ def make_type_global(input_deck): bank_source = particle_bank(0) bank_precursor = precursor_bank(0) - # Domain banks if needed - if input_deck.technique["domain_decomposition"]: - bank_domain_xp = particle_bank(input_deck.technique["domain_bank_size"]) - bank_domain_xn = particle_bank(input_deck.technique["domain_bank_size"]) - bank_domain_yp = particle_bank(input_deck.technique["domain_bank_size"]) - bank_domain_yn = particle_bank(input_deck.technique["domain_bank_size"]) - bank_domain_zp = particle_bank(input_deck.technique["domain_bank_size"]) - bank_domain_zn = particle_bank(input_deck.technique["domain_bank_size"]) - else: - bank_domain_xp = particle_bank(0) - bank_domain_xn = particle_bank(0) - bank_domain_yp = particle_bank(0) - bank_domain_yn = particle_bank(0) - bank_domain_zp = particle_bank(0) - bank_domain_zn = particle_bank(0) - bank_lost = particle_bank(0) # Particle tracker N_track = 0 if input_deck.setting["track_particle"]: @@ -1131,14 +1096,7 @@ def make_type_global(input_deck): ("bank_active", bank_active), ("bank_census", bank_census), ("bank_source", bank_source), - ("bank_domain_xp", bank_domain_xp), - ("bank_domain_xn", bank_domain_xn), - ("bank_domain_yp", bank_domain_yp), - ("bank_domain_yn", bank_domain_yn), - ("bank_domain_zp", bank_domain_zp), - ("bank_domain_zn", bank_domain_zn), ("bank_precursor", bank_precursor), - ("dd_idx", int64), ("k_eff", float64), ("k_cycle", float64, (N_cycle,)), ("k_avg", float64), diff --git a/test/regression/dd_slab_reed/answer.h5 b/test/regression/dd_slab_reed/answer.h5 deleted file mode 100644 index 8ea0bb48..00000000 Binary files a/test/regression/dd_slab_reed/answer.h5 and /dev/null differ diff --git a/test/regression/dd_slab_reed/input.py b/test/regression/dd_slab_reed/input.py deleted file mode 100644 index 167e102b..00000000 --- a/test/regression/dd_slab_reed/input.py +++ /dev/null @@ -1,51 +0,0 @@ -import numpy as np - -import mcdc - -# ============================================================================= -# Set model -# ============================================================================= -# Three slab layers with different materials -# Based on William H. Reed, NSE (1971), 46:2, 309-314, DOI: 10.13182/NSE46-309 - -# Set materials -m1 = mcdc.material(capture=np.array([50.0])) -m2 = mcdc.material(capture=np.array([5.0])) -m3 = mcdc.material(capture=np.array([0.0])) # Vacuum -m4 = mcdc.material(capture=np.array([0.1]), scatter=np.array([[0.9]])) - -# Set surfaces -s1 = mcdc.surface("plane-z", z=0.0, bc="reflective") -s2 = mcdc.surface("plane-z", z=2.0) -s3 = mcdc.surface("plane-z", z=3.0) -s4 = mcdc.surface("plane-z", z=5.0) -s5 = mcdc.surface("plane-z", z=8.0, bc="vacuum") - -# Set cells -mcdc.cell([+s1, -s2], m1) -mcdc.cell([+s2, -s3], m2) -mcdc.cell([+s3, -s4], m3) -mcdc.cell([+s4, -s5], m4) - -# ============================================================================= -# Set source -# ============================================================================= - -# Isotropic source in the absorbing medium -mcdc.source(z=[0.0, 2.0], isotropic=True, prob=50.0) - -# Isotropic source in the first half of the outermost medium, -# with 1/100 strength -mcdc.source(z=[5.0, 6.0], isotropic=True, prob=0.5) - -# ============================================================================= -# Set tally, setting, and run mcdc -# ============================================================================= - -mcdc.tally(scores=["flux"], z=np.linspace(0.0, 8.0, 81)) - -# Setting -mcdc.setting(N_particle=5000) -mcdc.domain_decomposition(z=np.linspace(0.0, 8.0, 5)) -# Run -mcdc.run() diff --git a/test/regression/run.py b/test/regression/run.py index 47bb181d..eac007ba 100644 --- a/test/regression/run.py +++ b/test/regression/run.py @@ -28,21 +28,11 @@ names = [item for item in os.listdir() if fnmatch.fnmatch(item, name)] names.sort() - # Remove skipped if specified if skip != "NONE": skips = [item for item in os.listdir() if fnmatch.fnmatch(item, skip)] for name in skips: names.remove(name) -# Skip cache if any -if "__pycache__" in names: - names.remove("__pycache__") -# Skip domain decomp tests unless there are 4 MPI processes -temp = names.copy() -for name in names: - if name[:3] == "dd_" and mpiexec != 4: - temp.remove(name) -names = temp # Data for each test printouts = [] @@ -54,6 +44,10 @@ # Run all tests for i, name in enumerate(names): + # Skip cache if any + if name == "__pycache__": + continue + print("\n[%i/%i] " % (i + 1, len(names)) + name) error_msgs.append([]) crashes.append(False)