Merge pull request #169 from shac170/good_dd

domain decomposition

.gitignore

@@ -51,3 +51,4 @@ pytestdebug.log
 docs/build
 docs/source/pythonapi/generated/
 
+*.csv

examples/fixed_source/slab_reed_dd/input.py

@@ -0,0 +1,51 @@
+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()

examples/fixed_source/slab_reed_dd/process.py

@@ -0,0 +1,123 @@
+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()

mcdc/__init__.py

@@ -23,6 +23,7 @@
     uq,
     print_card,
     reset_cards,
+    domain_decomposition,
 )
 from mcdc.main import run, prepare
 from mcdc.visualizer import visualize

mcdc/card.py

@@ -80,6 +80,12 @@ 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,

mcdc/constant.py

@@ -14,6 +14,7 @@
 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
@@ -38,6 +39,12 @@
 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)

mcdc/input_.py

@@ -1173,6 +1173,50 @@ 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,