Merge pull request #119 from spasmann/main

Parallelization of LDS and Additional Unit Tests

.gitignore

@@ -14,6 +14,7 @@ __pycache__
 *.swp
 
 # Output
+*output.h5
 output*.h5
 *output.h5
 *.png

mcdc/__init__.py

@@ -23,4 +23,4 @@
     print_card,
     reset_cards,
 )
-from mcdc.main import run
+from mcdc.main import run, prepare

mcdc/kernel.py

@@ -2602,12 +2602,10 @@ def prepare_qmc_particles(mcdc):
     QMC Low-Discrepency Sequence. Particles are added to the bank_source.
 
     """
-    # determine which portion of particles to loop through
+    # total number of particles
     N_particle = mcdc["setting"]["N_particle"]
+    # number of particles this processor will handle
     N_work = mcdc["mpi_work_size"]
-    rank = mcdc["mpi_rank"]
-    start = int(rank * N_work)
-    stop = int((rank + 1) * N_work)
 
     # low discrepency sequence
     lds = mcdc["technique"]["iqmc_lds"]
@@ -2627,7 +2625,7 @@ def prepare_qmc_particles(mcdc):
     za = mesh["z"][0]
     zb = mesh["z"][-1]
 
-    for n in range(start, stop):
+    for n in range(N_work):
         # Create new particle
         P_new = np.zeros(1, dtype=type_.particle_record)[0]
         P_new["rng_seed"] = 0

mcdc/main.py

@@ -363,35 +363,40 @@ def prepare():
         mcdc["technique"]["iqmc_mesh"]["z"] = input_deck.technique["iqmc_mesh"]["z"]
         mcdc["technique"]["iqmc_mesh"]["t"] = input_deck.technique["iqmc_mesh"]["t"]
         mcdc["technique"]["iqmc_generator"] = input_deck.technique["iqmc_generator"]
+        # variables to generate samples
+        scramble = mcdc["technique"]["iqmc_scramble"]
+        N_dim = mcdc["technique"]["iqmc_N_dim"]
+        seed = mcdc["technique"]["iqmc_seed"]
+        N = mcdc["setting"]["N_particle"]
+        size = MPI.COMM_WORLD.Get_size()
+        rank = MPI.COMM_WORLD.Get_rank()
+        N_work = math.ceil(N_particle / size)
+        # how many samples will we skip in the LDS
+        fast_forward = int((rank / size) * N)
         # generate lds
         if input_deck.technique["iqmc_generator"] == "sobol":
-            scramble = mcdc["technique"]["iqmc_scramble"]
-            N_dim = mcdc["technique"]["iqmc_N_dim"]
-            N = mcdc["setting"]["N_particle"]
             sampler = qmc.Sobol(d=N_dim, scramble=scramble)
-            m = math.ceil(math.log(N, 2))
-            mcdc["setting"]["N_particle"] = 2**m
-            mcdc["technique"]["iqmc_lds"] = sampler.random_base2(m=m)
-            # first row of an unscrambled Sobol sequence is all zeros
-            # and throws off some of the algorithms. So we add small amount
-            # to avoid this issue
-            # mcdc["technique"]["lds"][0,:] += 1e-6
-            mcdc["technique"]["iqmc_lds"][mcdc["technique"]["iqmc_lds"] == 0.0] += 1e-6
-            # lds is shape (2**m, d)
+            # skip first two entries in Sobol sequence because
+            # they map to x = 0.0 and ux = 0.0 respectively
+            sampler.fast_forward(2)
+            sampler.fast_forward(fast_forward)
+            mcdc["technique"]["iqmc_lds"] = sampler.random(N_work)
         if input_deck.technique["iqmc_generator"] == "halton":
-            scramble = mcdc["technique"]["iqmc_scramble"]
-            N_dim = mcdc["technique"]["iqmc_N_dim"]
-            seed = mcdc["technique"]["iqmc_seed"]
-            N = mcdc["setting"]["N_particle"]
             sampler = qmc.Halton(d=N_dim, scramble=scramble, seed=seed)
+            # skip first entry in Halton because it maps to x = 0.0
             sampler.fast_forward(1)
-            mcdc["technique"]["iqmc_lds"] = sampler.random(N)
+            sampler.fast_forward(fast_forward)
+            mcdc["technique"]["iqmc_lds"] = sampler.random(N_work)
         if input_deck.technique["iqmc_generator"] == "random":
-            seed = mcdc["technique"]["iqmc_seed"]
-            N_dim = mcdc["technique"]["iqmc_N_dim"]
-            N = mcdc["setting"]["N_particle"]
+            # this chunk of code uses the iqmc_seed to generate a number of
+            # seeds to be used  on each processor
+            # this way, each processor gets different samples, but if iQMC is run
+            # several times it will generate the same samples across runs
+            # 1e6 represents the maximum integer size generated
             np.random.seed(seed)
-            mcdc["technique"]["iqmc_lds"] = np.random.random((N, N_dim))
+            seeds = np.random.randint(1e6, size=size)
+            np.random.seed(seeds[rank])
+            mcdc["technique"]["iqmc_lds"] = np.random.random((N_work, N_dim))
 
     # =========================================================================
     # Derivative Source Method

mcdc/type_.py

@@ -550,7 +550,8 @@ def make_type_technique(N_particle, G, card):
     struct += [("iqmc_mesh", mesh)]
     # Low-discprenecy sequence
     # TODO: make N_dim an input setting
-    struct += [("iqmc_lds", float64, (N_particle, N_dim))]
+    N_work = math.ceil(N_particle / MPI.COMM_WORLD.Get_size())
+    struct += [("iqmc_lds", float64, (N_work, N_dim))]
 
     # Source
     struct += [("iqmc_source", float64, (Ng, Nt, Nx, Ny, Nz))]

test/unit/test_kernel.py

@@ -1,11 +1,74 @@
 import numpy as np
+import mcdc as MCDC
 from mcdc import type_
-from mcdc.kernel import rng
+from mcdc.main import closeout
+from mcdc.kernel import rng, AxV, FxV, HxV, preconditioner
 import mcdc.global_ as mcdc_
 
 input_deck = mcdc_.input_deck
 
 
+def iqmc_dummy_mcdc_variable():
+    """
+    This function returns the global "mcdc" container. Inputs are
+    taken from the kornreich eigenvalue problem.
+
+    """
+    # Set materials
+    m1 = MCDC.material(
+        capture=np.array([0.0]),
+        scatter=np.array([[0.9]]),
+        fission=np.array([0.1]),
+        nu_p=np.array([6.0]),
+    )
+    m2 = MCDC.material(
+        capture=np.array([0.68]),
+        scatter=np.array([[0.2]]),
+        fission=np.array([0.12]),
+        nu_p=np.array([2.5]),
+    )
+
+    # Set surfaces
+    s1 = MCDC.surface("plane-x", x=0.0, bc="vacuum")
+    s2 = MCDC.surface("plane-x", x=1.5)
+    s3 = MCDC.surface("plane-x", x=2.5, bc="vacuum")
+
+    # Set cells
+    MCDC.cell([+s1, -s2], m1)
+    MCDC.cell([+s2, -s3], m2)
+
+    # =============================================================================
+    # iQMC Parameters
+    # =============================================================================
+    N = 100
+    maxit = 10
+    tol = 1e-3
+    x = np.arange(0.0, 2.6, 0.1)
+    Nx = len(x) - 1
+    generator = "halton"
+    solver = "davidson"
+    fixed_source = np.zeros(Nx)
+    phi0 = np.ones((Nx))
+
+    # =============================================================================
+    # Set tally, setting, and run mcdc
+    # =============================================================================
+
+    MCDC.iQMC(
+        x=x,
+        fixed_source=fixed_source,
+        phi0=phi0,
+        maxitt=maxit,
+        tol=tol,
+        generator=generator,
+        eigenmode_solver=solver,
+    )
+    # Setting
+    MCDC.setting(N_particle=N)
+    MCDC.eigenmode()
+    return MCDC.prepare()
+
+
 def test_rn_basic():
     """
     Basic test routine for random number generator.
@@ -18,6 +81,8 @@ def test_rn_basic():
         Trans. Am. Nucl. Soc, 71, 202 (1994)
 
     """
+    MCDC.reset_cards()
+
     ref_data = np.array(
         (
             1,
@@ -28,36 +93,130 @@ def test_rn_basic():
         )
     )
 
-    # PRNG parameters from [1]
-    g = 2806196910506780709
-    c = 1
-    seed0 = 1
-    mod = 2**63
-
-    # Arbitrary constants for dummy mcdc container
-    G = 1
-    J = 1
-    Nmax_slice = 1
-    Nmax_surface = 1
-    Nmax_cell = 1
-    # need a psuedo material for mcdc container
-    input_deck.materials.append({"G": 1, "J": 1})
-
-    # Make types for dummy mcdc container
-    type_.make_type_material(G, J, 1)
-    type_.make_type_surface(Nmax_slice)
-    type_.make_type_cell(Nmax_surface)
-    type_.make_type_universe(Nmax_cell)
-    type_.make_type_lattice(input_deck.lattices)
-    type_.make_type_source(G)
-    type_.make_type_tally(1, input_deck.tally)
-    type_.make_type_technique(0, 1, input_deck.technique)
-    type_.make_type_global(input_deck)
-
-    # The dummy container
-    mcdc = np.zeros(1, dtype=type_.global_)[0]
+    mcdc = iqmc_dummy_mcdc_variable()
 
     # run through the first five seeds (1-5)
     for i in range(5):
-        rng(mcdc["setting"])
         assert mcdc["setting"]["rng_seed"] == ref_data[i]
+        rng(mcdc["setting"])
+
+
+def test_AxV_linearity():
+    """
+    AxV is the linear operator used for GMRES in iQMC.
+
+    Linear operators must satisfy conditions of additivity and multiplicity
+    defined as:
+            - Additivity: f(x+y) = f(x) + f(y)
+            - Multiplicity: f(cx) = cf(x)
+
+    We can test both properties with:
+            - f(a*x + b*y) = a*f(x) + b*f(y)
+    """
+    MCDC.reset_cards()
+
+    mcdc = iqmc_dummy_mcdc_variable()
+
+    size = mcdc["technique"]["iqmc_flux"].size
+    np.random.seed(123456)
+    a = np.random.random()
+    b = np.random.random()
+    x = np.random.random((size,))
+    y = np.random.random((size,))
+    rhs = np.zeros((size,))
+
+    F1 = AxV((a * x + b * y), rhs, mcdc)
+    F2 = a * AxV(x, rhs, mcdc) + b * AxV(y, rhs, mcdc)
+    assert np.allclose(F1, F2, rtol=1e-10)
+
+
+def test_FxV_linearity():
+    """
+    FxV is a linear operator used for Davidson in iQMC.
+
+    Linear operators must satisfy conditions of additivity and multiplicity
+    defined as:
+            - Additivity: f(x+y) = f(x) + f(y)
+            - Multiplicity: f(cx) = cf(x)
+
+    We can test both properties with:
+            - f(a*x + b*y) = a*f(x) + b*f(y)
+    """
+    MCDC.reset_cards()
+
+    mcdc = iqmc_dummy_mcdc_variable()
+
+    size = mcdc["technique"]["iqmc_flux"].size
+    np.random.seed(123456)
+    a = np.random.random()
+    b = np.random.random()
+    x = np.random.random((size, 1))
+    y = np.random.random((size, 1))
+
+    F1 = FxV((a * x + b * y), mcdc)
+    F2 = a * FxV(x, mcdc) + b * FxV(y, mcdc)
+    assert np.allclose(F1, F2, rtol=1e-10)
+
+
+def test_HxV_linearity():
+    """
+    HxV is a linear operator used for Davidson in iQMC.
+
+    Linear operators must satisfy conditions of additivity and multiplicity
+    defined as:
+            - Additivity: f(x+y) = f(x) + f(y)
+            - Multiplicity: f(cx) = cf(x)
+
+    We can test both properties with:
+            - f(a*x + b*y) = a*f(x) + b*f(y)
+    """
+    MCDC.reset_cards()
+
+    mcdc = iqmc_dummy_mcdc_variable()
+
+    size = mcdc["technique"]["iqmc_flux"].size
+    np.random.seed(123456)
+    a = np.random.random()
+    b = np.random.random()
+    x = np.random.random((size, 1))
+    y = np.random.random((size, 1))
+
+    F1 = HxV((a * x + b * y), mcdc)
+    F2 = a * HxV(x, mcdc) + b * HxV(y, mcdc)
+    assert np.allclose(F1, F2, rtol=1e-10)
+
+
+def test_preconditioner_linearity():
+    """
+    perconditioner is a linear operator used for Davidson in iQMC.
+
+    Linear operators must satisfy conditions of additivity and multiplicity
+    defined as:
+            - Additivity: f(x+y) = f(x) + f(y)
+            - Multiplicity: f(cx) = cf(x)
+
+    We can test both properties with:
+            - f(a*x + b*y) = a*f(x) + b*f(y)
+    """
+    MCDC.reset_cards()
+
+    mcdc = iqmc_dummy_mcdc_variable()
+
+    size = mcdc["technique"]["iqmc_flux"].size
+    np.random.seed(123456)
+    a = np.random.random()
+    b = np.random.random()
+    x = np.random.random((size, 1))
+    y = np.random.random((size, 1))
+
+    F1 = preconditioner((a * x + b * y), mcdc)
+    F2 = a * preconditioner(x, mcdc) + b * preconditioner(y, mcdc)
+    assert np.allclose(F1, F2, rtol=1e-10)
+
+
+# if __name__ == "__main__":
+# test_rn_basic()
+# test_AxV_linearity()
+#     test_FxV_linearity()
+#     test_HxV_linearity()
+#     test_preconditioner_linearity()