This page is based on a Jupyter notebook.

import os… (click to toggle)

import os
from pathlib import Path

import numpy as np
import ogstools as ot
import pyvista as pv
from IPython.display import Markdown, display
from numpy.random import default_rng
from scipy.spatial import cKDTree

… (click to toggle)

out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
    out_dir.mkdir(parents=True)
orig_dir = Path.cwd()

Reactive transport and retention in stochastically generated DFNs: A PorePy–OpenGeoSys Workflow

Protecting the environment over geological timescales requires isolating hazardous substances (e.g., radionuclides and industrial chemicals) from the biosphere. Accurate prediction of solute transport is essential for the safe design of deep geological repositories and for long-term groundwater contamination protection.

Why is solute retention critical?

Natural retention processes — matrix diffusion and sorption — serve as critical barriers that delay contaminant migration, helping to prevent their spread through the environment.

Matrix Diffusion: The movement of solutes from regions of high to low concentration within a solid matrix, governed by concentration gradients and diffusion coefficients.
Sorption: The process by which solutes are retained on or in a solid, encompassing adsorption (solutes attaching to the surface) and absorption (solutes entering the solid matrix). This process is influenced by surface chemistry, solute properties, and environmental conditions, and is quantified by the distribution coefficient.

How to study the solute retention?

Tracer tests are a primary method for investigating subsurface solute transport and retention. In these experiments, a known tracer is injected into the formation, and its concentration is monitored over time at a recovery location.

Different tracers are used to isolate transport mechanisms:

Nonsorbing tracers (e.g., tritiated water, bromide) reflect advective transport and matrix diffusion, providing information on physical transport mechanisms without chemical retention effects.
Sorption-active tracers (e.g., strontium, cesium) additionally capture chemical retention through sorption.

How to model the transport in fractured rock?

In crystalline rock formations, groundwater flow occurs primarily through discrete fracture networks (DFN), as the rock matrix itself exhibits low permeability. DFN models simulate flow and solute migration by treating fractures as discrete hydraulic features.

How to numerically model it?

To investigate solute transport and retention in fractured rocks, a DFN is generated using PorePy with stochastic fracture placement algorithms. The resulting DFN geometry is processed through OGSTools to create OpenGeoSys (OGS)-compatible input files, enabling the setup and execution of a Hydro-Component (HC) simulation for flow and solute transport, and retention in realistic fractured media.

Problem description

alt text

OpenGeosys

See DFNbyPorePy.py how to generate the DFN mesh and geometry data required by this notebook.

Preparing the mesh (`.vtu`) for OpenGeosys simulation

We generate s clean, standardized .vtu mesh containing only MaterialIDs as cell data — fully tailored for OpenGeoSys simulations in fractured media.

Loads the .vtu file into a mesh object for processing.
If MaterialIDs are missing, they’re generated from subdomain_id, starting at 0 — OGS uses them to identify physical regions.
Removes all other cell and point data, keeping only what’s necessary for OGS.
Saves the meshes ready for use in an OGS .prj setup.

DFN_2D = ot.mesh.read("mixed_dimensional_grid_2.vtu")… (click to toggle)

DFN_2D = ot.mesh.read("mixed_dimensional_grid_2.vtu")
if "MaterialIDs" not in DFN_2D.cell_data:
    DFN_2D["MaterialIDs"] = (
        DFN_2D["subdomain_id"] - DFN_2D["subdomain_id"].min()
    ).astype(np.int32)

for key in list(DFN_2D.cell_data):
    if key != "MaterialIDs":
        DFN_2D.cell_data.remove(key)

for key in list(DFN_2D.point_data):
    DFN_2D.point_data.remove(key)


DFN_2D.save(f"{out_dir}/mixed_dimensional_grid_2.vtu")

os.chdir(out_dir)… (click to toggle)

os.chdir(out_dir)

ot.cli().reviseMesh(i="mixed_dimensional_grid_2.vtu", o="temp.vtu")
ot.cli().NodeReordering(i="temp.vtu", o="mixed_dimensional_grid_2.vtu")
ot.cli().checkMesh("mixed_dimensional_grid_2.vtu", v=True)

info: Mesh read: 10340 nodes, 20662 elements.
info: Simplifying the mesh...
warning: Property MaterialIDs exists but does not have the requested mesh item type node.
warning: Property MaterialIDs exists but does not have the requested type float.
warning: Property MaterialIDs exists but does not have the requested type double.
info: Revised mesh: 10340 nodes, 20662 elements.
info: Reordering nodes... 
info: Method: Reversing order of nodes unless it is considered correct by the OGS6 standard, i.e. such that det(J) > 0, where J is the Jacobian of the global-to-local coordinate transformation.

info: Corrected 0 elements.
info: VTU file written.
info: Memory size: 3 MiB
info: Time for reading: 0.0274005 s
info: Axis aligned bounding box: 	x [-1.77636e-15, 100) (extent 100)
	y [0, 100) (extent 100)
	z [0, 100) (extent 100)
info: Min/max edge lengths: [0.0336252, 6.85976]
info: Number of elements in the mesh:
info: 	Triangles: 20662
info: Mesh Quality Control:
info: Looking for unused nodes...
info: Found 0 potentially collapsible nodes.
info: Testing mesh element geometry:
info: 10360 elements found with wrong node order.
info: No elements found with zero volume.

info: No elements found with non coplanar nodes.

info: No elements found with non-convex geometry.

info: 10360 elements found with wrong node order.
ElementIDs: 0, 2, 5, 6, 8, 9, 13, 14, 15, 17, 18, 20, 22, 24, 26, 27, 28, 30, 33, 35, 37, 39, 41, 42, 43, 45, 47, 50, 51, 53, 57, 59, 62, 63, 66, 67, 70, 72, 73, 77, 80, 83, 85, 86, 88, 92, 95, 97, 98, 100, 101, 104, 105, 107, 109, 110, 113, 115, 117, 119, 120, 124, 128, 132, 138, 142, 146, 151, 152, 154, 155, 156, 158, 159, 160, 161, 162, 164, 166, 167, 168, 170, 171, 178, 181, 185, 188, 193, 196, 198, 200, 202, 203, 204, 206, 208, 212, 214, 216, 220, 223, 226, 230, 232, 235, 237, 245, 249, 253, 257, 261, 267, 269, 273, 276, 277, 278, 279, 281, 282, 283, 284, 286, 287, 288, 291, 292, 293, 295, 296, 297, 299, 301, 303, 304, 307, 309, 311, 313, 315, 317, 319, 321, 322, 325, 326, 328, 330, 331, 333, 335, 337, 339, 343, 346, 347, 349, 351, 354, 355, 359, 360, 363, 366, 367, 369, 370, 372, 373, 375, 377, 380, 384, 386, 387, 389, 391, 392, 395, 398, 401, 402, 405, 406, 409, 410, 412, 413, 417, 420, 424, 428, 431, 433, 434, 435, 436, 437, 439, 440, 441, 444, 446, 449, 451, 454, 456, 459, 461, 464, 465, 467, 469, 471, 472, 473, 475, 478, 480, 481, 483, 486, 488, 490, 493, 495, 497, 501, 505, 508, 510, 511, 512, 513, 515, 516, 517, 519, 520, 521, 523, 525, 527, 528, 530, 532, 534, 536, 539, 541, 542, 544, 547, 548, 550, 552, 553, 556, 558, 559, 560, 562, 564, 566, 568, 570, 571, 572, 573, 576, 580, 582, 585, 586, 588, 589, 592, 593, 596, 597, 598, 599, 601, 603, 608, 609, 612, 613, 615, 616, 617, 620, 621, 625, 628, 629, 631, 634, 635, 637, 638, 639, 641, 642, 643, 647, 649, 651, 653, 655, 657, 658, 659, 660, 662, 664, 666, 668, 671, 673, 675, 676, 678, 680, 681, 682, 683, 685, 686, 688, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699, 702, 703, 705, 713, 715, 716, 719, 721, 722, 727, 728, 730, 731, 732, 736, 737, 738, 741, 742, 743, 747, 748, 749, 750, 753, 757, 758, 760, 761, 762, 764, 769, 770, 771, 772, 775, 778, 779, 784, 785, 786, 787, 789, 790, 791, 796, 797, 798, 800, 803, 805, 806, 807, 808, 812, 813, 818, 819, 820, 822, 823, 829, 830, 831, 832, 836, 841, 843, 844, 845, 846, 847, 849, 852, 854, 856, 857, 858, 861, 863, 865, 866, 867, 868, 871, 872, 878, 879, 880, 885, 886, 887, 891, 894, 895, 896, 897, 898, 899, 901, 902, 906, 907, 908, 912, 914, 915, 916, 917, 919, 920, 926, 930, 934, 935, 937, 939, 943, 944, 946, 947, 948, 950, 952, 956, 959, 960, 961, 963, 964, 965, 968, 969, 970, 972, 973, 974, 977, 978, 979, 980, 985, 986, 987, 990, 991, 992, 997, 998, 999, 1002, 1003, 1004, 1009, 1010, 1011, 1014, 1015, 1016, 1017, 1021, 1022, 1027, 1028, 1032, 1033, 1034, 1035, 1037, 1038, 1043, 1044, 1045, 1046, 1048, 1049, 1050, 1051, 1053, 1054, 1058, 1059, 1061, 1062, 1063, 1065, 1067, 1071, 1072, 1073, 1074, 1075, 1080, 1081, 1086, 1087, 1088, 1091, 1092, 1093, 1098, 1099, 1103, 1104, 1105, 1106, 1110, 1111, 1115, 1116, 1117, 1118, 1122, 1123, 1124, 1127, 1133, 1134, 1135, 1136, 1141, 1143, 1145, 1147, 1148, 1149, 1150, 1151, 1152, 1153, 1156, 1160, 1163, 1165, 1166, 1169, 1170, 1171, 1175, 1176, 1177, 1181, 1182, 1183, 1186, 1187, 1189, 1191, 1193, 1196, 1197, 1198, 1202, 1203, 1204, 1205, 1206, 1208, 1209, 1210, 1215, 1216, 1217, 1218, 1219, 1220, 1221, 1226, 1227, 1228, 1233, 1235, 1236, 1237, 1238, 1239, 1240, 1245, 1248, 1256, 1258, 1259, 1260, 1261, 1262, 1263, 1264, 1267, 1277, 1278, 1279, 1280, 1281, 1283, 1286, 1287, 1288, 1291, 1292, 1293, 1297, 1298, 1299, 1300, 1301, 1303, 1304, 1305, 1306, 1307, 1308, 1312, 1313, 1317, 1318, 1319, 1320, 1327, 1334, 1335, 1339, 1340, 1342, 1343, 1344, 1350, 1351, 1353, 1354, 1356, 1357, 1359, 1362, 1364, 1367, 1371, 1375, 1376, 1377, 1378, 1383, 1384, 1386, 1387, 1391, 1392, 1393, 1394, 1395, 1397, 1398, 1403, 1404, 1405, 1409, 1410, 1411, 1412, 1413, 1415, 1416, 1420, 1424, 1425, 1426, 1427, 1428, 1430, 1436, 1437, 1438, 1441, 1442, 1443, 1448, 1449, 1450, 1451, 1453, 1454, 1461, 1462, 1463, 1466, 1467, 1468, 1473, 1474, 1475, 1477, 1478, 1482, 1483, 1484, 1487, 1488, 1489, 1490, 1494, 1495, 1496, 1497, 1501, 1504, 1506, 1511, 1512, 1515, 1516, 1517, 1518, 1521, 1522, 1527, 1528, 1529, 1530, 1535, 1537, 1538, 1540, 1541, 1545, 1547, 1548, 1554, 1555, 1558, 1559, 1560, 1561, 1562, 1563, 1565, 1566, 1570, 1571, 1572, 1577, 1578, 1582, 1583, 1584, 1585, 1588, 1589, 1590, 1595, 1597, 1598, 1603, 1604, 1607, 1608, 1611, 1612, 1613, 1617, 1618, 1619, 1622, 1623, 1624, 1625, 1629, 1630, 1635, 1636, 1637, 1638, 1641, 1645, 1647, 1648, 1649, 1650, 1652, 1657, 1658, 1660, 1661, 1662, 1667, 1669, 1670, 1672, 1673, 1674, 1678, 1679, 1680, 1681, 1682, 1685, 1687, 1692, 1694, 1695, 1696, 1697, 1703, 1704, 1705, 1709, 1710, 1715, 1716, 1717, 1718, 1720, 1725, 1726, 1727, 1728, 1729, 1730, 1736, 1737, 1738, 1739, 1743, 1744, 1745, 1746, 1748, 1749, 1750, 1755, 1756, 1757, 1758, 1764, 1765, 1766, 1767, 1771, 1772, 1777, 1778, 1779, 1782, 1786, 1787, 1788, 1789, 1790, 1791, 1793, 1794, 1796, 1801, 1802, 1803, 1805, 1808, 1809, 1813, 1814, 1816, 1819, 1821, 1823, 1824, 1828, 1829, 1833, 1834, 1835, 1836, 1841, 1842, 1843, 1845, 1848, 1852, 1854, 1856, 1859, 1860, 1863, 1865, 1869, 1871, 1874, 1877, 1878, 1880, 1881, 1882, 1887, 1888, 1891, 1894, 1899, 1900, 1901, 1903, 1904, 1905, 1906, 1909, 1910, 1913, 1916, 1918, 1921, 1924, 1926, 1928, 1929, 1934, 1935, 1937, 1939, 1941, 1944, 1947, 1950, 1951, 1952, 1953, 1955, 1960, 1961, 1965, 1966, 1968, 1969, 1972, 1974, 1976, 1978, 1979, 1982, 1984, 1985, 1987, 1990, 1992, 1993, 1997, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2020, 2022, 2025, 2026, 2027, 2028, 2029, 2030, 2032, 2033, 2035, 2036, 2037, 2039, 2040, 2043, 2045, 2046, 2049, 2050, 2052, 2059, 2060, 2061, 2062, 2064, 2065, 2071, 2072, 2079, 2080, 2082, 2085, 2090, 2092, 2094, 2095, 2100, 2102, 2106, 2110, 2112, 2114, 2116, 2118, 2120, 2122, 2125, 2128, 2130, 2133, 2135, 2136, 2143, 2145, 2147, 2148, 2149, 2150, 2152, 2155, 2157, 2158, 2159, 2162, 2164, 2166, 2169, 2171, 2173, 2174, 2176, 2177, 2180, 2181, 2182, 2183, 2185, 2190, 2191, 2193, 2194, 2197, 2198, 2199, 2202, 2203, 2204, 2206, 2207, 2208, 2212, 2215, 2216, 2217, 2220, 2221, 2223, 2227, 2228, 2229, 2231, 2232, 2233, 2234, 2235, 2239, 2240, 2241, 2242, 2243, 2246, 2249, 2250, 2251, 2252, 2254, 2258, 2264, 2267, 2268, 2270, 2272, 2273, 2274, 2275, 2276, 2277, 2279, 2284, 2286, 2287, 2288, 2289, 2290, 2294, 2295, 2296, 2298, 2300, 2301, 2305, 2306, 2307, 2308, 2311, 2315, 2317, 2319, 2320, 2324, 2325, 2326, 2327, 2329, 2332, 2334, 2337, 2338, 2339, 2340, 2343, 2344, 2347, 2349, 2352, 2356, 2359, 2360, 2361, 2363, 2364, 2367, 2368, 2371, 2373, 2375, 2379, 2380, 2381, 2384, 2388, 2389, 2393, 2394, 2396, 2397, 2399, 2402, 2405, 2406, 2408, 2409, 2413, 2414, 2417, 2418, 2420, 2422, 2425, 2427, 2428, 2429, 2433, 2436, 2437, 2438, 2440, 2443, 2445, 2446, 2447, 2449, 2451, 2453, 2457, 2458, 2462, 2468, 2472, 2474, 2475, 2476, 2480, 2483, 2485, 2486, 2487, 2489, 2491, 2493, 2495, 2497, 2498, 2500, 2502, 2503, 2505, 2506, 2508, 2509, 2511, 2513, 2514, 2516, 2517, 2524, 2528, 2530, 2534, 2535, 2541, 2545, 2547, 2551, 2553, 2557, 2561, 2562, 2563, 2565, 2567, 2568, 2569, 2570, 2571, 2573, 2574, 2576, 2580, 2582, 2586, 2588, 2592, 2594, 2598, 2604, 2608, 2613, 2616, 2618, 2619, 2620, 2621, 2622, 2624, 2625, 2626, 2628, 2629, 2633, 2636, 2638, 2639, 2644, 2645, 2648, 2649, 2650, 2651, 2653, 2654, 2656, 2659, 2661, 2662, 2663, 2666, 2667, 2668, 2670, 2672, 2674, 2675, 2676, 2678, 2681, 2682, 2685, 2687, 2688, 2689, 2690, 2692, 2695, 2697, 2698, 2700, 2702, 2703, 2705, 2707, 2709, 2710, 2711, 2714, 2716, 2718, 2720, 2722, 2723, 2724, 2726, 2727, 2728, 2732, 2734, 2736, 2738, 2739, 2741, 2743, 2745, 2747, 2749, 2751, 2753, 2755, 2756, 2757, 2759, 2763, 2767, 2769, 2773, 2779, 2781, 2783, 2785, 2787, 2788, 2790, 2793, 2796, 2800, 2802, 2804, 2805, 2806, 2807, 2809, 2811, 2813, 2814, 2815, 2817, 2819, 2820, 2824, 2826, 2827, 2828, 2830, 2831, 2833, 2837, 2839, 2840, 2841, 2843, 2844, 2845, 2848, 2850, 2854, 2855, 2857, 2859, 2861, 2862, 2865, 2867, 2868, 2871, 2873, 2875, 2877, 2878, 2879, 2881, 2882, 2884, 2885, 2886, 2887, 2889, 2890, 2892, 2894, 2895, 2896, 2898, 2900, 2903, 2906, 2908, 2912, 2914, 2917, 2919, 2921, 2923, 2924, 2925, 2928, 2930, 2934, 2935, 2938, 2940, 2941, 2943, 2944, 2945, 2948, 2949, 2952, 2955, 2959, 2960, 2962, 2966, 2968, 2969, 2971, 2975, 2978, 2979, 2980, 2982, 2984, 2987, 2989, 2991, 2993, 2995, 2998, 3004, 3007, 3009, 3010, 3011, 3016, 3017, 3020, 3021, 3023, 3026, 3028, 3029, 3030, 3032, 3034, 3035, 3039, 3041, 3042, 3045, 3050, 3053, 3054, 3058, 3061, 3067, 3069, 3070, 3071, 3072, 3079, 3081, 3083, 3084, 3085, 3088, 3093, 3094, 3095, 3096, 3102, 3103, 3104, 3107, 3111, 3112, 3113, 3114, 3115, 3120, 3121, 3122, 3123, 3126, 3127, 3130, 3131, 3136, 3141, 3142, 3143, 3144, 3145, 3146, 3148, 3149, 3154, 3155, 3162, 3164, 3165, 3166, 3167, 3172, 3173, 3174, 3175, 3176, 3180, 3181, 3182, 3183, 3184, 3185, 3186, 3187, 3192, 3193, 3197, 3198, 3199, 3205, 3211, 3213, 3214, 3217, 3218, 3219, 3223, 3224, 3226, 3227, 3232, 3233, 3234, 3237, 3240, 3241, 3242, 3243, 3247, 3250, 3251, 3252, 3254, 3259, 3261, 3263, 3264, 3265, 3270, 3272, 3273, 3274, 3275, 3276, 3277, 3280, 3284, 3285, 3286, 3287, 3288, 3293, 3294, 3296, 3297, 3299, 3306, 3307, 3308, 3309, 3310, 3311, 3316, 3317, 3318, 3320, 3323, 3325, 3329, 3330, 3331, 3334, 3335, 3336, 3341, 3342, 3344, 3346, 3351, 3355, 3356, 3358, 3360, 3361, 3364, 3366, 3367, 3368, 3371, 3374, 3375, 3377, 3378, 3379, 3380, 3382, 3383, 3385, 3386, 3389, 3390, 3391, 3392, 3393, 3396, 3397, 3398, 3400, 3402, 3403, 3404, 3406, 3412, 3413, 3416, 3417, 3418, 3419, 3423, 3426, 3427, 3429, 3431, 3432, 3433, 3435, 3436, 3438, 3439, 3440, 3442, 3445, 3447, 3448, 3450, 3451, 3452, 3458, 3459, 3460, 3461, 3465, 3466, 3469, 3470, 3471, 3478, 3479, 3480, 3481, 3483, 3484, 3485, 3486, 3492, 3493, 3494, 3495, 3499, 3500, 3501, 3505, 3506, 3507, 3508, 3511, 3512, 3513, 3514, 3515, 3517, 3518, 3520, 3522, 3523, 3524, 3526, 3527, 3530, 3532, 3534, 3535, 3537, 3540, 3544, 3545, 3546, 3548, 3550, 3554, 3555, 3557, 3558, 3560, 3563, 3564, 3566, 3567, 3569, 3570, 3573, 3574, 3577, 3580, 3585, 3588, 3590, 3591, 3594, 3595, 3597, 3599, 3600, 3601, 3602, 3603, 3611, 3614, 3617, 3618, 3621, 3622, 3628, 3629, 3630, 3631, 3632, 3633, 3634, 3638, 3639, 3644, 3645, 3646, 3649, 3650, 3651, 3656, 3657, 3658, 3660, 3661, 3665, 3666, 3668, 3669, 3670, 3673, 3674, 3676, 3681, 3682, 3683, 3684, 3685, 3686, 3687, 3692, 3693, 3694, 3697, 3699, 3701, 3702, 3705, 3709, 3710, 3711, 3716, 3717, 3721, 3722, 3723, 3724, 3725, 3726, 3727, 3728, 3733, 3734, 3735, 3739, 3741, 3743, 3745, 3746, 3747, 3748, 3750, 3751, 3755, 3756, 3757, 3760, 3761, 3762, 3767, 3768, 3769, 3771, 3773, 3774, 3779, 3780, 3781, 3785, 3786, 3787, 3788, 3789, 3792, 3793, 3798, 3799, 3805, 3806, 3808, 3818, 3819, 3820, 3822, 3824, 3825, 3826, 3827, 3829, 3830, 3831, 3832, 3833, 3834, 3835, 3836, 3837, 3838, 3839, 3842, 3843, 3844, 3850, 3852, 3855, 3858, 3859, 3860, 3862, 3863, 3866, 3869, 3873, 3877, 3878, 3880, 3881, 3886, 3887, 3895, 3896, 3901, 3903, 3904, 3907, 3908, 3913, 3914, 3917, 3918, 3920, 3921, 3923, 3924, 3928, 3930, 3931, 3934, 3935, 3938, 3942, 3945, 3946, 3948, 3951, 3953, 3955, 3956, 3957, 3958, 3963, 3964, 3965, 3969, 3970, 3971, 3972, 3975, 3977, 3978, 3979, 3980, 3981, 3982, 3986, 3987, 3988, 3989, 3992, 3994, 3995, 3997, 3998, 3999, 4002, 4003, 4007, 4010, 4011, 4013, 4014, 4017, 4021, 4022, 4026, 4027, 4028, 4035, 4036, 4038, 4040, 4042, 4046, 4047, 4048, 4049, 4051, 4053, 4055, 4057, 4061, 4063, 4067, 4068, 4069, 4071, 4076, 4077, 4079, 4081, 4083, 4085, 4086, 4088, 4089, 4090, 4091, 4092, 4093, 4095, 4096, 4097, 4098, 4101, 4103, 4105, 4106, 4108, 4110, 4111, 4112, 4117, 4118, 4119, 4120, 4122, 4125, 4128, 4129, 4130, 4132, 4136, 4138, 4140, 4142, 4143, 4144, 4145, 4148, 4149, 4151, 4154, 4155, 4158, 4159, 4162, 4163, 4165, 4168, 4170, 4172, 4173, 4174, 4176, 4178, 4182, 4183, 4184, 4185, 4187, 4191, 4195, 4199, 4201, 4202, 4204, 4205, 4207, 4209, 4210, 4212, 4214, 4216, 4218, 4219, 4221, 4225, 4227, 4228, 4232, 4234, 4235, 4237, 4238, 4239, 4240, 4242, 4243, 4245, 4247, 4250, 4251, 4253, 4257, 4259, 4260, 4262, 4263, 4268, 4270, 4272, 4274, 4276, 4279, 4280, 4282, 4284, 4286, 4288, 4294, 4298, 4302, 4306, 4309, 4313, 4317, 4319, 4320, 4321, 4322, 4323, 4324, 4326, 4327, 4328, 4330, 4332, 4338, 4340, 4341, 4345, 4346, 4347, 4349, 4350, 4352, 4353, 4354, 4355, 4356, 4357, 4359, 4360, 4361, 4362, 4365, 4368, 4369, 4370, 4371, 4372, 4374, 4376, 4380, 4383, 4390, 4394, 4396, 4400, 4401, 4402, 4406, 4409, 4410, 4412, 4415, 4417, 4419, 4420, 4422, 4424, 4425, 4426, 4428, 4429, 4430, 4432, 4435, 4436, 4437, 4439, 4442, 4444, 4447, 4450, 4452, 4454, 4455, 4456, 4458, 4459, 4461, 4463, 4464, 4465, 4466, 4468, 4469, 4470, 4472, 4474, 4476, 4477, 4479, 4480, 4483, 4485, 4488, 4490, 4491, 4493, 4494, 4495, 4498, 4505, 4506, 4507, 4509, 4512, 4513, 4517, 4518, 4520, 4522, 4525, 4526, 4528, 4529, 4531, 4533, 4534, 4536, 4537, 4540, 4541, 4542, 4543, 4545, 4548, 4551, 4552, 4553, 4555, 4556, 4558, 4559, 4561, 4563, 4571, 4572, 4577, 4578, 4580, 4581, 4584, 4585, 4586, 4588, 4589, 4591, 4594, 4596, 4598, 4603, 4604, 4606, 4608, 4610, 4612, 4615, 4617, 4619, 4622, 4624, 4626, 4629, 4631, 4633, 4635, 4637, 4640, 4642, 4644, 4646, 4649, 4651, 4654, 4656, 4658, 4660, 4661, 4663, 4665, 4668, 4669, 4671, 4673, 4676, 4677, 4682, 4684, 4685, 4687, 4689, 4692, 4694, 4698, 4700, 4701, 4703, 4706, 4708, 4709, 4711, 4714, 4716, 4718, 4719, 4720, 4722, 4723, 4725, 4727, 4728, 4730, 4732, 4734, 4738, 4740, 4742, 4746, 4754, 4755, 4758, 4759, 4762, 4764, 4767, 4769, 4770, 4771, 4773, 4775, 4777, 4778, 4781, 4784, 4785, 4787, 4788, 4792, 4795, 4798, 4800, 4802, 4805, 4807, 4808, 4810, 4813, 4815, 4818, 4819, 4822, 4823, 4825, 4827, 4829, 4832, 4833, 4835, 4837, 4840, 4843, 4845, 4846, 4847, 4848, 4851, 4857, 4859, 4860, 4861, 4863, 4864, 4865, 4867, 4871, 4872, 4873, 4876, 4878, 4879, 4880, 4882, 4883, 4884, 4886, 4887, 4891, 4896, 4898, 4901, 4902, 4905, 4906, 4908, 4909, 4910, 4911, 4913, 4914, 4917, 4921, 4922, 4924, 4926, 4928, 4929, 4930, 4932, 4933, 4934, 4936, 4938, 4941, 4942, 4945, 4946, 4947, 4951, 4952, 4953, 4955, 4956, 4959, 4960, 4961, 4962, 4965, 4969, 4970, 4972, 4974, 4976, 4978, 4983, 4985, 4988, 4991, 4993, 4996, 4997, 5000, 5005, 5007, 5008, 5009, 5012, 5016, 5018, 5020, 5021, 5023, 5027, 5028, 5029, 5030, 5037, 5038, 5039, 5040, 5041, 5045, 5046, 5047, 5048, 5052, 5053, 5057, 5058, 5059, 5060, 5064, 5065, 5069, 5070, 5071, 5072, 5076, 5077, 5078, 5080, 5081, 5087, 5088, 5089, 5090, 5091, 5092, 5096, 5097, 5098, 5099, 5105, 5111, 5112, 5113, 5117, 5118, 5119, 5120, 5121, 5122, 5123, 5124, 5125, 5126, 5131, 5132, 5136, 5137, 5138, 5139, 5142, 5143, 5144, 5149, 5150, 5151, 5153, 5154, 5158, 5161, 5162, 5163, 5165, 5166, 5170, 5171, 5172, 5175, 5176, 5177, 5178, 5181, 5182, 5183, 5188, 5189, 5190, 5191, 5194, 5198, 5199, 5200, 5201, 5202, 5205, 5206, 5208, 5214, 5216, 5217, 5218, 5219, 5220, 5224, 5225, 5226, 5227, 5233, 5234, 5235, 5239, 5240, 5243, 5244, 5245, 5251, 5252, 5255, 5256, 5257, 5258, 5259, 5261, 5262, 5263, 5267, 5268, 5269, 5273, 5278, 5282, 5283, 5284, 5285, 5286, 5288, 5289, 5290, 5294, 5295, 5296, 5300, 5301, 5302, 5303, 5304, 5305, 5306, 5312, 5319, 5320, 5321, 5322, 5324, 5328, 5331, 5332, 5335, 5336, 5337, 5341, 5344, 5345, 5346, 5351, 5353, 5355, 5356, 5357, 5360, 5362, 5363, 5364, 5367, 5368, 5371, 5372, 5375, 5377, 5378, 5379, 5384, 5385, 5386, 5387, 5388, 5389, 5395, 5396, 5397, 5398, 5401, 5404, 5405, 5406, 5407, 5411, 5412, 5413, 5416, 5418, 5420, 5421, 5423, 5425, 5427, 5430, 5432, 5434, 5435, 5436, 5438, 5439, 5440, 5442, 5446, 5448, 5458, 5459, 5460, 5464, 5465, 5466, 5467, 5471, 5472, 5476, 5477, 5478, 5479, 5480, 5486, 5488, 5491, 5495, 5496, 5497, 5498, 5503, 5504, 5505, 5507, 5508, 5509, 5521, 5522, 5523, 5524, 5528, 5531, 5532, 5533, 5534, 5536, 5540, 5541, 5542, 5546, 5547, 5548, 5549, 5550, 5551, 5552, 5559, 5563, 5565, 5566, 5567, 5571, 5574, 5575, 5576, 5580, 5581, 5582, 5589, 5590, 5591, 5592, 5593, 5594, 5595, 5596, 5597, 5598, 5599, 5600, 5601, 5607, 5608, 5609, 5613, 5616, 5617, 5620, 5621, 5622, 5624, 5626, 5627, 5632, 5635, 5636, 5637, 5639, 5642, 5643, 5644, 5646, 5648, 5651, 5652, 5656, 5657, 5658, 5659, 5661, 5663, 5664, 5667, 5668, 5669, 5671, 5672, 5675, 5680, 5681, 5683, 5684, 5685, 5690, 5691, 5692, 5695, 5698, 5699, 5700, 5704, 5705, 5709, 5711, 5712, 5714, 5715, 5716, 5717, 5718, 5719, 5721, 5722, 5725, 5726, 5727, 5728, 5729, 5731, 5732, 5733, 5736, 5737, 5740, 5741, 5742, 5744, 5745, 5746, 5749, 5750, 5751, 5753, 5755, 5756, 5758, 5759, 5763, 5765, 5766, 5768, 5771, 5773, 5775, 5777, 5778, 5780, 5781, 5783, 5786, 5787, 5789, 5790, 5794, 5795, 5796, 5798, 5801, 5806, 5807, 5810, 5811, 5812, 5813, 5815, 5818, 5820, 5824, 5830, 5832, 5837, 5838, 5839, 5840, 5844, 5845, 5847, 5848, 5850, 5856, 5857, 5858, 5859, 5862, 5865, 5867, 5868, 5870, 5872, 5874, 5875, 5877, 5878, 5883, 5884, 5885, 5888, 5889, 5890, 5895, 5896, 5897, 5898, 5899, 5904, 5905, 5906, 5909, 5913, 5914, 5915, 5916, 5917, 5922, 5923, 5924, 5925, 5926, 5930, 5931, 5932, 5935, 5940, 5941, 5942, 5946, 5947, 5948, 5949, 5952, 5956, 5958, 5960, 5965, 5966, 5967, 5968, 5969, 5971, 5973, 5974, 5975, 5978, 5980, 5985, 5988, 5989, 5993, 5997, 6000, 6001, 6003, 6005, 6006, 6007, 6010, 6016, 6018, 6019, 6021, 6022, 6023, 6024, 6025, 6029, 6031, 6032, 6037, 6039, 6040, 6042, 6043, 6045, 6047, 6049, 6050, 6053, 6054, 6059, 6062, 6064, 6066, 6067, 6068, 6069, 6071, 6072, 6073, 6075, 6079, 6082, 6085, 6086, 6087, 6090, 6093, 6094, 6095, 6096, 6099, 6101, 6102, 6103, 6104, 6105, 6108, 6110, 6111, 6112, 6113, 6117, 6118, 6119, 6120, 6124, 6125, 6126, 6127, 6132, 6133, 6137, 6139, 6140, 6141, 6142, 6143, 6144, 6149, 6150, 6151, 6152, 6154, 6155, 6157, 6158, 6160, 6164, 6165, 6166, 6168, 6170, 6171, 6172, 6173, 6176, 6177, 6179, 6183, 6186, 6190, 6192, 6193, 6194, 6195, 6196, 6198, 6199, 6205, 6206, 6207, 6209, 6211, 6212, 6213, 6216, 6218, 6220, 6224, 6225, 6226, 6229, 6230, 6231, 6236, 6237, 6239, 6240, 6241, 6243, 6244, 6248, 6250, 6251, 6252, 6253, 6255, 6257, 6261, 6262, 6263, 6264, 6268, 6271, 6274, 6276, 6278, 6280, 6281, 6282, 6283, 6286, 6288, 6289, 6296, 6300, 6303, 6304, 6305, 6307, 6315, 6319, 6320, 6321, 6322, 6325, 6326, 6327, 6328, 6330, 6333, 6335, 6337, 6339, 6340, 6341, 6344, 6347, 6348, 6349, 6352, 6353, 6354, 6355, 6362, 6363, 6365, 6366, 6372, 6374, 6375, 6377, 6380, 6382, 6383, 6385, 6386, 6387, 6389, 6390, 6391, 6395, 6397, 6398, 6399, 6401, 6403, 6405, 6407, 6409, 6413, 6415, 6416, 6417, 6418, 6419, 6421, 6422, 6423, 6424, 6426, 6427, 6428, 6432, 6436, 6437, 6438, 6439, 6440, 6442, 6443, 6447, 6448, 6454, 6458, 6460, 6461, 6464, 6468, 6469, 6470, 6471, 6474, 6475, 6477, 6479, 6480, 6482, 6484, 6485, 6488, 6491, 6493, 6494, 6495, 6497, 6498, 6500, 6502, 6504, 6506, 6509, 6510, 6512, 6513, 6514, 6515, 6517, 6518, 6520, 6522, 6523, 6524, 6526, 6528, 6529, 6530, 6531, 6533, 6534, 6538, 6540, 6542, 6544, 6546, 6547, 6551, 6553, 6555, 6558, 6560, 6562, 6563, 6565, 6568, 6570, 6576, 6580, 6581, 6582, 6586, 6588, 6592, 6596, 6600, 6604, 6613, 6617, 6621, 6627, 6628, 6634, 6636, 6640, 6642, 6643, 6644, 6646, 6647, 6648, 6650, 6651, 6652, 6654, 6655, 6659, 6661, 6663, 6667, 6671, 6675, 6679, 6685, 6686, 6687, 6688, 6689, 6691, 6692, 6694, 6696, 6698, 6701, 6702, 6706, 6709, 6711, 6715, 6716, 6718, 6720, 6721, 6722, 6725, 6727, 6729, 6730, 6732, 6733, 6734, 6736, 6740, 6744, 6746, 6748, 6751, 6753, 6756, 6758, 6760, 6763, 6767, 6768, 6772, 6773, 6774, 6776, 6777, 6780, 6781, 6784, 6785, 6787, 6790, 6792, 6793, 6795, 6797, 6800, 6801, 6802, 6804, 6805, 6807, 6808, 6813, 6816, 6819, 6821, 6822, 6824, 6825, 6826, 6828, 6829, 6831, 6833, 6835, 6836, 6838, 6840, 6841, 6843, 6845, 6847, 6848, 6850, 6853, 6856, 6860, 6864, 6866, 6868, 6870, 6873, 6875, 6878, 6880, 6881, 6883, 6884, 6885, 6890, 6892, 6893, 6894, 6896, 6897, 6898, 6901, 6903, 6905, 6907, 6908, 6910, 6913, 6916, 6918, 6921, 6923, 6924, 6927, 6928, 6931, 6932, 6935, 6938, 6942, 6944, 6945, 6947, 6948, 6950, 6953, 6955, 6958, 6960, 6962, 6963, 6965, 6967, 6968, 6971, 6973, 6974, 6976, 6979, 6983, 6986, 6987, 6989, 6990, 6992, 6994, 6995, 6996, 6998, 7002, 7004, 7006, 7008, 7009, 7010, 7013, 7014, 7017, 7020, 7021, 7024, 7026, 7028, 7030, 7033, 7036, 7037, 7039, 7041, 7042, 7044, 7047, 7048, 7050, 7051, 7053, 7054, 7055, 7056, 7057, 7060, 7061, 7062, 7064, 7067, 7068, 7069, 7071, 7072, 7074, 7076, 7077, 7078, 7079, 7080, 7081, 7083, 7084, 7085, 7086, 7093, 7096, 7097, 7098, 7099, 7100, 7105, 7107, 7109, 7111, 7112, 7113, 7114, 7118, 7119, 7120, 7121, 7122, 7123, 7126, 7127, 7128, 7133, 7134, 7135, 7136, 7138, 7145, 7146, 7147, 7149, 7150, 7154, 7155, 7156, 7157, 7160, 7162, 7163, 7165, 7166, 7167, 7170, 7171, 7172, 7175, 7179, 7181, 7183, 7184, 7185, 7187, 7188, 7189, 7191, 7192, 7193, 7194, 7196, 7201, 7203, 7207, 7209, 7210, 7216, 7217, 7218, 7222, 7223, 7225, 7228, 7229, 7230, 7231, 7235, 7238, 7239, 7240, 7243, 7247, 7248, 7249, 7250, 7253, 7254, 7256, 7257, 7262, 7263, 7264, 7266, 7268, 7269, 7271, 7272, 7275, 7276, 7280, 7281, 7282, 7283, 7286, 7287, 7288, 7289, 7291, 7294, 7295, 7299, 7300, 7301, 7302, 7308, 7309, 7310, 7313, 7315, 7317, 7318, 7319, 7323, 7326, 7328, 7330, 7331, 7332, 7333, 7336, 7340, 7346, 7347, 7352, 7353, 7357, 7361, 7366, 7367, 7369, 7370, 7371, 7372, 7373, 7377, 7378, 7379, 7382, 7383, 7384, 7385, 7388, 7389, 7390, 7391, 7395, 7396, 7399, 7401, 7402, 7403, 7404, 7406, 7407, 7410, 7414, 7415, 7416, 7417, 7424, 7426, 7428, 7433, 7434, 7437, 7439, 7440, 7444, 7447, 7449, 7450, 7451, 7453, 7457, 7458, 7462, 7463, 7464, 7465, 7469, 7472, 7473, 7474, 7476, 7481, 7482, 7484, 7485, 7490, 7491, 7492, 7494, 7495, 7499, 7500, 7501, 7502, 7504, 7511, 7512, 7513, 7514, 7515, 7516, 7517, 7518, 7519, 7523, 7524, 7526, 7527, 7528, 7529, 7530, 7531, 7532, 7533, 7540, 7542, 7543, 7545, 7551, 7552, 7554, 7555, 7557, 7559, 7562, 7568, 7570, 7577, 7578, 7580, 7583, 7584, 7588, 7589, 7591, 7592, 7595, 7600, 7602, 7603, 7604, 7606, 7607, 7608, 7609, 7610, 7612, 7614, 7616, 7617, 7618, 7619, 7620, 7622, 7623, 7624, 7625, 7626, 7630, 7631, 7632, 7634, 7636, 7639, 7640, 7642, 7644, 7646, 7648, 7650, 7652, 7654, 7655, 7659, 7660, 7662, 7663, 7664, 7665, 7666, 7670, 7673, 7675, 7681, 7683, 7684, 7689, 7690, 7693, 7694, 7695, 7697, 7702, 7703, 7705, 7706, 7708, 7710, 7713, 7714, 7716, 7717, 7720, 7722, 7723, 7725, 7726, 7727, 7729, 7730, 7735, 7736, 7740, 7741, 7742, 7743, 7744, 7748, 7751, 7752, 7755, 7756, 7757, 7763, 7764, 7765, 7766, 7770, 7771, 7774, 7776, 7777, 7778, 7779, 7780, 7785, 7788, 7789, 7790, 7793, 7795, 7797, 7798, 7800, 7801, 7803, 7804, 7806, 7807, 7811, 7814, 7815, 7817, 7819, 7821, 7823, 7828, 7829, 7830, 7833, 7837, 7838, 7840, 7841, 7842, 7844, 7845, 7846, 7849, 7851, 7854, 7856, 7858, 7859, 7860, 7861, 7863, 7866, 7867, 7868, 7869, 7870, 7873, 7875, 7879, 7882, 7885, 7886, 7888, 7891, 7896, 7897, 7899, 7901, 7902, 7904, 7906, 7907, 7909, 7910, 7911, 7913, 7915, 7918, 7920, 7922, 7923, 7925, 7926, 7928, 7930, 7932, 7933, 7936, 7937, 7940, 7941, 7944, 7945, 7947, 7949, 7950, 7951, 7954, 7955, 7957, 7961, 7962, 7965, 7967, 7968, 7970, 7971, 7972, 7975, 7977, 7979, 7981, 7983, 7985, 7986, 7988, 7989, 7991, 7993, 7996, 7998, 8000, 8002, 8003, 8005, 8006, 8007, 8009, 8012, 8014, 8017, 8019, 8021, 8022, 8023, 8025, 8027, 8030, 8032, 8033, 8035, 8038, 8042, 8043, 8044, 8046, 8048, 8049, 8050, 8052, 8053, 8054, 8056, 8057, 8058, 8060, 8061, 8062, 8064, 8065, 8066, 8069, 8070, 8072, 8073, 8074, 8076, 8077, 8078, 8080, 8081, 8082, 8083, 8086, 8088, 8091, 8093, 8094, 8098, 8099, 8100, 8101, 8102, 8104, 8105, 8106, 8107, 8109, 8110, 8111, 8112, 8114, 8116, 8117, 8119, 8120, 8123, 8125, 8127, 8129, 8131, 8133, 8135, 8137, 8139, 8141, 8142, 8143, 8145, 8147, 8150, 8153, 8154, 8158, 8161, 8162, 8165, 8167, 8168, 8170, 8171, 8175, 8176, 8177, 8181, 8183, 8184, 8185, 8187, 8190, 8191, 8195, 8196, 8197, 8199, 8204, 8206, 8208, 8210, 8211, 8212, 8214, 8218, 8219, 8221, 8224, 8226, 8229, 8232, 8233, 8234, 8236, 8238, 8243, 8244, 8247, 8250, 8251, 8255, 8258, 8260, 8262, 8263, 8265, 8268, 8270, 8271, 8273, 8274, 8276, 8277, 8281, 8283, 8284, 8286, 8287, 8288, 8290, 8291, 8292, 8294, 8297, 8300, 8302, 8305, 8307, 8309, 8311, 8313, 8316, 8318, 8319, 8320, 8322, 8324, 8326, 8329, 8330, 8333, 8336, 8338, 8340, 8342, 8344, 8346, 8350, 8352, 8354, 8356, 8358, 8359, 8361, 8363, 8364, 8366, 8368, 8370, 8371, 8374, 8375, 8376, 8377, 8379, 8380, 8381, 8383, 8385, 8388, 8389, 8392, 8394, 8396, 8398, 8399, 8402, 8406, 8407, 8408, 8409, 8411, 8412, 8413, 8415, 8416, 8417, 8420, 8422, 8424, 8425, 8427, 8428, 8429, 8432, 8433, 8436, 8437, 8439, 8440, 8442, 8444, 8446, 8450, 8454, 8456, 8458, 8462, 8464, 8466, 8470, 8472, 8474, 8475, 8477, 8480, 8482, 8483, 8484, 8486, 8490, 8492, 8493, 8495, 8496, 8497, 8498, 8500, 8503, 8505, 8508, 8512, 8513, 8515, 8517, 8519, 8522, 8524, 8527, 8529, 8530, 8531, 8533, 8535, 8537, 8540, 8541, 8542, 8545, 8547, 8550, 8553, 8556, 8559, 8561, 8563, 8565, 8566, 8569, 8571, 8575, 8578, 8580, 8583, 8585, 8587, 8588, 8590, 8592, 8593, 8595, 8597, 8600, 8602, 8604, 8607, 8609, 8611, 8614, 8615, 8616, 8617, 8618, 8619, 8621, 8622, 8623, 8625, 8626, 8627, 8629, 8630, 8631, 8634, 8638, 8642, 8644, 8646, 8650, 8652, 8654, 8659, 8662, 8664, 8666, 8671, 8673, 8674, 8676, 8679, 8684, 8687, 8689, 8692, 8693, 8696, 8698, 8701, 8705, 8707, 8710, 8713, 8715, 8718, 8720, 8722, 8723, 8726, 8727, 8728, 8730, 8731, 8734, 8735, 8737, 8739, 8740, 8742, 8743, 8745, 8747, 8749, 8751, 8754, 8755, 8756, 8759, 8761, 8763, 8765, 8767, 8769, 8771, 8773, 8774, 8775, 8777, 8779, 8782, 8784, 8787, 8789, 8791, 8792, 8794, 8797, 8800, 8802, 8804, 8806, 8810, 8811, 8812, 8813, 8815, 8816, 8817, 8818, 8820, 8821, 8822, 8824, 8827, 8829, 8830, 8833, 8834, 8836, 8838, 8839, 8842, 8844, 8845, 8846, 8849, 8851, 8852, 8853, 8854, 8857, 8858, 8859, 8862, 8868, 8872, 8874, 8876, 8879, 8881, 8884, 8886, 8887, 8888, 8889, 8893, 8894, 8895, 8901, 8902, 8903, 8904, 8910, 8912, 8913, 8914, 8915, 8916, 8918, 8919, 8920, 8921, 8924, 8925, 8926, 8927, 8929, 8932, 8936, 8938, 8939, 8941, 8942, 8946, 8949, 8950, 8951, 8955, 8956, 8960, 8961, 8962, 8963, 8967, 8970, 8971, 8972, 8973, 8974, 8975, 8979, 8980, 8982, 8983, 8986, 8988, 8989, 8990, 8991, 8995, 8996, 8999, 9004, 9005, 9006, 9007, 9008, 9009, 9013, 9015, 9017, 9018, 9019, 9020, 9021, 9023, 9025, 9030, 9031, 9033, 9035, 9036, 9044, 9045, 9048, 9049, 9050, 9051, 9058, 9059, 9064, 9068, 9069, 9070, 9071, 9073, 9077, 9078, 9080, 9081, 9085, 9086, 9088, 9091, 9096, 9100, 9101, 9102, 9104, 9106, 9108, 9109, 9111, 9112, 9113, 9115, 9116, 9118, 9119, 9120, 9123, 9127, 9129, 9132, 9133, 9134, 9135, 9141, 9142, 9143, 9144, 9145, 9149, 9150, 9154, 9155, 9156, 9157, 9163, 9164, 9165, 9166, 9170, 9171, 9174, 9175, 9176, 9177, 9182, 9183, 9184, 9185, 9188, 9191, 9193, 9198, 9200, 9201, 9203, 9204, 9205, 9208, 9214, 9215, 9216, 9217, 9218, 9221, 9222, 9223, 9224, 9229, 9230, 9232, 9233, 9236, 9237, 9238, 9239, 9240, 9241, 9242, 9243, 9244, 9249, 9250, 9251, 9258, 9260, 9263, 9267, 9268, 9269, 9270, 9271, 9273, 9278, 9282, 9283, 9284, 9285, 9290, 9291, 9292, 9294, 9295, 9298, 9300, 9304, 9307, 9310, 9312, 9313, 9314, 9316, 9317, 9319, 9320, 9321, 9322, 9323, 9324, 9326, 9328, 9330, 9331, 9335, 9337, 9338, 9339, 9341, 9343, 9344, 9347, 9348, 9349, 9351, 9353, 9356, 9357, 9358, 9359, 9360, 9361, 9364, 9367, 9370, 9372, 9373, 9378, 9380, 9384, 9385, 9386, 9387, 9389, 9393, 9395, 9396, 9399, 9403, 9405, 9406, 9407, 9408, 9414, 9415, 9416, 9417, 9421, 9422, 9423, 9424, 9425, 9429, 9430, 9431, 9432, 9436, 9437, 9441, 9442, 9443, 9444, 9448, 9449, 9450, 9452, 9453, 9457, 9459, 9460, 9461, 9462, 9465, 9468, 9469, 9470, 9471, 9475, 9476, 9480, 9481, 9482, 9483, 9487, 9490, 9491, 9492, 9495, 9496, 9497, 9502, 9503, 9504, 9507, 9508, 9509, 9510, 9514, 9515, 9519, 9520, 9521, 9522, 9529, 9531, 9532, 9533, 9534, 9540, 9541, 9542, 9547, 9550, 9551, 9553, 9554, 9555, 9559, 9562, 9565, 9566, 9567, 9573, 9574, 9575, 9576, 9580, 9581, 9582, 9583, 9584, 9589, 9590, 9592, 9593, 9594, 9597, 9599, 9601, 9602, 9603, 9605, 9607, 9608, 9609, 9611, 9612, 9619, 9620, 9623, 9624, 9628, 9629, 9630, 9635, 9636, 9643, 9644, 9645, 9648, 9649, 9650, 9651, 9655, 9656, 9661, 9662, 9663, 9666, 9667, 9668, 9673, 9674, 9675, 9678, 9679, 9680, 9685, 9686, 9687, 9691, 9693, 9696, 9697, 9699, 9700, 9701, 9702, 9706, 9709, 9710, 9714, 9715, 9716, 9717, 9721, 9722, 9725, 9726, 9727, 9729, 9730, 9733, 9735, 9736, 9737, 9739, 9740, 9741, 9742, 9745, 9747, 9749, 9750, 9751, 9752, 9753, 9754, 9755, 9760, 9761, 9763, 9764, 9767, 9768, 9769, 9770, 9774, 9775, 9779, 9780, 9781, 9782, 9786, 9788, 9792, 9795, 9799, 9800, 9801, 9805, 9808, 9809, 9810, 9811, 9812, 9813, 9816, 9817, 9818, 9819, 9824, 9825, 9827, 9828, 9829, 9830, 9835, 9836, 9837, 9838, 9841, 9845, 9846, 9847, 9848, 9850, 9853, 9857, 9858, 9859, 9860, 9866, 9870, 9871, 9876, 9877, 9878, 9880, 9881, 9887, 9888, 9889, 9890, 9891, 9893, 9894, 9896, 9897, 9898, 9899, 9902, 9903, 9904, 9909, 9910, 9911, 9913, 9917, 9918, 9919, 9920, 9922, 9927, 9928, 9929, 9930, 9931, 9937, 9939, 9941, 9942, 9943, 9948, 9952, 9953, 9957, 9958, 9959, 9961, 9962, 9966, 9971, 9974, 9975, 9977, 9978, 9979, 9981, 9986, 9987, 9989, 9995, 9997, 9999, 10001, 10002, 10003, 10004, 10005, 10006, 10019, 10022, 10023, 10024, 10028, 10029, 10033, 10034, 10035, 10040, 10041, 10042, 10043, 10045, 10052, 10053, 10054, 10057, 10058, 10060, 10065, 10067, 10068, 10072, 10075, 10078, 10080, 10081, 10082, 10083, 10084, 10085, 10090, 10091, 10092, 10094, 10095, 10099, 10100, 10102, 10103, 10104, 10107, 10111, 10113, 10114, 10115, 10116, 10120, 10121, 10125, 10126, 10127, 10128, 10132, 10133, 10137, 10138, 10139, 10140, 10146, 10147, 10148, 10149, 10153, 10156, 10157, 10158, 10160, 10161, 10164, 10168, 10169, 10171, 10172, 10173, 10179, 10180, 10181, 10182, 10183, 10184, 10185, 10189, 10191, 10192, 10193, 10194, 10195, 10196, 10198, 10199, 10200, 10202, 10204, 10207, 10208, 10209, 10210, 10211, 10212, 10213, 10214, 10215, 10222, 10225, 10226, 10227, 10231, 10234, 10235, 10236, 10237, 10238, 10239, 10243, 10244, 10245, 10246, 10250, 10251, 10252, 10255, 10256, 10260, 10262, 10263, 10264, 10265, 10266, 10270, 10274, 10277, 10279, 10281, 10286, 10287, 10288, 10289, 10291, 10292, 10293, 10295, 10296, 10297, 10299, 10300, 10303, 10304, 10305, 10312, 10314, 10315, 10318, 10319, 10320, 10321, 10322, 10324, 10331, 10332, 10333, 10339, 10340, 10341, 10342, 10343, 10346, 10347, 10348, 10352, 10353, 10354, 10358, 10360, 10367, 10368, 10370, 10371, 10372, 10373, 10374, 10378, 10379, 10380, 10385, 10386, 10387, 10388, 10389, 10390, 10397, 10398, 10399, 10400, 10401, 10402, 10403, 10404, 10405, 10406, 10407, 10412, 10415, 10416, 10417, 10423, 10432, 10436, 10437, 10439, 10440, 10442, 10447, 10449, 10452, 10453, 10454, 10456, 10459, 10460, 10465, 10466, 10467, 10469, 10470, 10471, 10472, 10473, 10476, 10477, 10478, 10479, 10480, 10481, 10486, 10487, 10492, 10493, 10494, 10497, 10498, 10499, 10500, 10501, 10502, 10503, 10510, 10513, 10515, 10516, 10517, 10518, 10520, 10522, 10525, 10527, 10529, 10532, 10533, 10538, 10539, 10540, 10544, 10545, 10546, 10550, 10551, 10552, 10555, 10556, 10557, 10558, 10559, 10561, 10563, 10565, 10566, 10570, 10572, 10573, 10576, 10577, 10580, 10582, 10585, 10586, 10587, 10588, 10592, 10595, 10596, 10598, 10600, 10604, 10607, 10608, 10610, 10611, 10612, 10613, 10614, 10615, 10618, 10619, 10620, 10627, 10628, 10633, 10634, 10639, 10640, 10641, 10642, 10643, 10646, 10647, 10648, 10650, 10652, 10656, 10657, 10661, 10662, 10663, 10664, 10668, 10670, 10672, 10674, 10676, 10680, 10683, 10687, 10688, 10689, 10690, 10693, 10695, 10697, 10699, 10700, 10701, 10704, 10705, 10708, 10711, 10712, 10715, 10716, 10718, 10720, 10722, 10723, 10724, 10726, 10727, 10732, 10733, 10735, 10736, 10737, 10739, 10740, 10741, 10746, 10747, 10748, 10751, 10752, 10755, 10757, 10764, 10767, 10768, 10770, 10771, 10774, 10776, 10777, 10779, 10781, 10784, 10787, 10790, 10791, 10795, 10799, 10802, 10804, 10806, 10808, 10809, 10811, 10817, 10818, 10820, 10822, 10824, 10826, 10827, 10828, 10830, 10835, 10838, 10839, 10840, 10843, 10845, 10847, 10848, 10849, 10850, 10856, 10857, 10862, 10865, 10868, 10869, 10870, 10872, 10873, 10875, 10876, 10877, 10879, 10880, 10882, 10883, 10884, 10888, 10889, 10890, 10891, 10892, 10895, 10897, 10898, 10902, 10907, 10908, 10909, 10913, 10914, 10919, 10920, 10925, 10928, 10931, 10936, 10940, 10942, 10943, 10945, 10946, 10947, 10948, 10950, 10951, 10952, 10955, 10959, 10960, 10961, 10962, 10964, 10965, 10968, 10969, 10972, 10973, 10976, 10977, 10978, 10980, 10982, 10983, 10984, 10986, 10987, 10988, 10990, 10991, 10993, 10994, 10995, 10998, 11002, 11004, 11007, 11012, 11014, 11018, 11026, 11027, 11030, 11034, 11035, 11036, 11040, 11041, 11042, 11045, 11046, 11047, 11051, 11052, 11053, 11054, 11061, 11062, 11064, 11067, 11068, 11070, 11073, 11074, 11080, 11086, 11087, 11088, 11089, 11092, 11095, 11096, 11097, 11098, 11099, 11101, 11103, 11105, 11106, 11107, 11108, 11109, 11110, 11112, 11113, 11116, 11117, 11118, 11119, 11120, 11121, 11122, 11123, 11124, 11126, 11127, 11131, 11132, 11133, 11134, 11135, 11138, 11139, 11140, 11141, 11143, 11144, 11146, 11150, 11154, 11155, 11156, 11158, 11160, 11161, 11162, 11165, 11168, 11169, 11170, 11171, 11172, 11175, 11176, 11180, 11181, 11182, 11183, 11187, 11188, 11190, 11191, 11192, 11194, 11195, 11197, 11198, 11199, 11201, 11202, 11206, 11207, 11208, 11209, 11210, 11214, 11216, 11217, 11219, 11221, 11223, 11225, 11228, 11229, 11230, 11234, 11235, 11236, 11241, 11244, 11246, 11247, 11248, 11250, 11251, 11253, 11254, 11256, 11257, 11258, 11259, 11263, 11265, 11266, 11268, 11273, 11276, 11283, 11284, 11285, 11294, 11295, 11296, 11300, 11301, 11302, 11303, 11304, 11306, 11308, 11312, 11313, 11314, 11316, 11318, 11320, 11321, 11323, 11325, 11326, 11331, 11332, 11333, 11334, 11336, 11337, 11338, 11340, 11343, 11345, 11348, 11350, 11351, 11353, 11354, 11357, 11359, 11360, 11361, 11362, 11364, 11366, 11369, 11372, 11373, 11376, 11377, 11378, 11380, 11381, 11382, 11383, 11384, 11386, 11389, 11390, 11393, 11395, 11396, 11397, 11400, 11403, 11405, 11407, 11410, 11412, 11414, 11415, 11417, 11418, 11421, 11422, 11425, 11426, 11428, 11429, 11431, 11433, 11434, 11437, 11438, 11441, 11442, 11443, 11445, 11448, 11449, 11451, 11454, 11456, 11457, 11458, 11460, 11461, 11463, 11464, 11465, 11466, 11467, 11468, 11470, 11471, 11472, 11474, 11475, 11476, 11477, 11478, 11480, 11481, 11482, 11483, 11484, 11486, 11487, 11488, 11494, 11495, 11499, 11506, 11508, 11511, 11512, 11513, 11515, 11516, 11517, 11519, 11520, 11521, 11523, 11524, 11526, 11527, 11528, 11530, 11531, 11532, 11534, 11535, 11536, 11538, 11539, 11541, 11542, 11543, 11546, 11547, 11548, 11549, 11550, 11552, 11553, 11556, 11558, 11560, 11565, 11569, 11572, 11573, 11577, 11579, 11581, 11583, 11585, 11587, 11589, 11595, 11596, 11597, 11598, 11600, 11601, 11604, 11607, 11608, 11609, 11611, 11612, 11615, 11617, 11618, 11621, 11622, 11624, 11625, 11628, 11629, 11630, 11631, 11633, 11634, 11637, 11640, 11642, 11643, 11646, 11647, 11648, 11649, 11651, 11653, 11655, 11657, 11658, 11660, 11661, 11665, 11666, 11668, 11670, 11672, 11673, 11676, 11677, 11680, 11686, 11687, 11689, 11690, 11693, 11694, 11696, 11697, 11699, 11701, 11705, 11707, 11708, 11709, 11710, 11711, 11714, 11715, 11716, 11718, 11719, 11721, 11730, 11733, 11740, 11743, 11745, 11750, 11751, 11753, 11754, 11755, 11759, 11761, 11763, 11764, 11765, 11766, 11767, 11768, 11769, 11771, 11772, 11777, 11778, 11782, 11783, 11784, 11785, 11791, 11792, 11793, 11798, 11799, 11800, 11802, 11803, 11807, 11808, 11809, 11813, 11814, 11815, 11816, 11817, 11818, 11821, 11822, 11823, 11824, 11827, 11828, 11829, 11834, 11835, 11836, 11839, 11840, 11841, 11846, 11847, 11848, 11850, 11851, 11859, 11860, 11861, 11862, 11863, 11864, 11865, 11866, 11871, 11872, 11873, 11875, 11876, 11882, 11883, 11884, 11885, 11889, 11890, 11892, 11893, 11898, 11902, 11903, 11904, 11909, 11910, 11911, 11912, 11913, 11916, 11917, 11921, 11922, 11923, 11926, 11927, 11928, 11933, 11934, 11935, 11940, 11942, 11943, 11944, 11947, 11949, 11950, 11955, 11958, 11959, 11960, 11961, 11962, 11963, 11965, 11966, 11971, 11972, 11974, 11976, 11981, 11984, 11985, 11989, 11991, 11992, 11997, 11998, 11999, 12003, 12005, 12008, 12011, 12012, 12013, 12014, 12015, 12016, 12021, 12023, 12024, 12025, 12026, 12029, 12032, 12033, 12034, 12035, 12040, 12041, 12042, 12047, 12048, 12049, 12050, 12053, 12054, 12055, 12059, 12060, 12064, 12065, 12066, 12067, 12072, 12074, 12075, 12076, 12077, 12078, 12079, 12084, 12086, 12087, 12088, 12089, 12092, 12093, 12095, 12096, 12097, 12099, 12105, 12106, 12107, 12112, 12113, 12114, 12115, 12117, 12118, 12120, 12123, 12128, 12129, 12132, 12134, 12135, 12136, 12138, 12141, 12142, 12144, 12146, 12147, 12151, 12153, 12154, 12159, 12163, 12164, 12165, 12166, 12168, 12170, 12172, 12175, 12176, 12178, 12180, 12181, 12183, 12184, 12187, 12191, 12193, 12194, 12199, 12201, 12203, 12208, 12209, 12211, 12213, 12217, 12219, 12220, 12221, 12222, 12224, 12227, 12229, 12230, 12231, 12233, 12235, 12236, 12238, 12239, 12241, 12243, 12246, 12247, 12251, 12252, 12255, 12256, 12257, 12258, 12259, 12260, 12267, 12270, 12271, 12272, 12273, 12274, 12277, 12278, 12281, 12284, 12286, 12287, 12288, 12290, 12292, 12296, 12298, 12300, 12302, 12303, 12305, 12306, 12307, 12310, 12311, 12312, 12313, 12316, 12317, 12318, 12321, 12322, 12325, 12326, 12328, 12329, 12331, 12334, 12335, 12337, 12339, 12341, 12342, 12344, 12346, 12348, 12350, 12353, 12355, 12356, 12359, 12360, 12361, 12367, 12368, 12371, 12372, 12373, 12375, 12378, 12381, 12382, 12383, 12385, 12387, 12390, 12392, 12395, 12398, 12399, 12401, 12405, 12408, 12410, 12411, 12412, 12414, 12415, 12416, 12418, 12419, 12420, 12422, 12423, 12424, 12426, 12427, 12428, 12430, 12431, 12432, 12434, 12435, 12436, 12439, 12441, 12442, 12443, 12444, 12445, 12447, 12448, 12449, 12451, 12452, 12454, 12456, 12457, 12464, 12468, 12471, 12472, 12473, 12475, 12476, 12477, 12478, 12480, 12482, 12483, 12484, 12486, 12487, 12489, 12490, 12491, 12492, 12493, 12495, 12498, 12499, 12500, 12502, 12503, 12505, 12507, 12509, 12510, 12511, 12512, 12514, 12515, 12516, 12517, 12518, 12520, 12521, 12522, 12523, 12526, 12530, 12533, 12536, 12539, 12541, 12542, 12543, 12547, 12550, 12553, 12557, 12560, 12562, 12565, 12567, 12568, 12570, 12571, 12573, 12577, 12578, 12580, 12582, 12584, 12586, 12589, 12590, 12592, 12594, 12596, 12597, 12599, 12601, 12603, 12605, 12606, 12607, 12609, 12610, 12613, 12617, 12618, 12619, 12620, 12621, 12622, 12624, 12627, 12630, 12631, 12635, 12637, 12641, 12642, 12644, 12646, 12648, 12649, 12651, 12653, 12654, 12655, 12656, 12658, 12659, 12661, 12662, 12663, 12664, 12665, 12666, 12670, 12673, 12675, 12679, 12681, 12682, 12684, 12686, 12688, 12689, 12690, 12692, 12693, 12694, 12696, 12697, 12699, 12701, 12702, 12704, 12708, 12709, 12710, 12712, 12713, 12715, 12718, 12720, 12722, 12724, 12728, 12729, 12731, 12733, 12735, 12737, 12741, 12744, 12745, 12748, 12749, 12751, 12754, 12756, 12758, 12759, 12762, 12763, 12765, 12767, 12770, 12772, 12773, 12775, 12777, 12780, 12782, 12784, 12786, 12787, 12790, 12791, 12793, 12795, 12796, 12797, 12802, 12806, 12810, 12812, 12814, 12817, 12819, 12821, 12822, 12824, 12825, 12828, 12831, 12834, 12836, 12838, 12839, 12841, 12843, 12846, 12847, 12852, 12856, 12858, 12860, 12862, 12863, 12865, 12867, 12870, 12871, 12873, 12875, 12880, 12884, 12886, 12888, 12891, 12894, 12895, 12896, 12900, 12904, 12906, 12907, 12909, 12913, 12915, 12916, 12918, 12921, 12924, 12925, 12928, 12930, 12931, 12933, 12934, 12939, 12943, 12947, 12949, 12951, 12952, 12954, 12956, 12957, 12958, 12960, 12963, 12965, 12967, 12970, 12971, 12973, 12974, 12976, 12977, 12981, 12983, 12985, 12986, 12987, 12989, 12990, 12993, 12995, 12996, 12998, 13001, 13003, 13004, 13005, 13007, 13009, 13011, 13012, 13013, 13016, 13017, 13020, 13021, 13023, 13025, 13026, 13027, 13029, 13031, 13032, 13034, 13038, 13039, 13042, 13043, 13045, 13047, 13048, 13050, 13051, 13052, 13053, 13056, 13059, 13062, 13063, 13067, 13069, 13070, 13072, 13075, 13077, 13078, 13080, 13081, 13082, 13083, 13084, 13086, 13087, 13089, 13090, 13094, 13096, 13098, 13100, 13101, 13102, 13103, 13104, 13106, 13109, 13112, 13113, 13114, 13117, 13118, 13120, 13121, 13123, 13125, 13128, 13129, 13130, 13132, 13134, 13135, 13138, 13139, 13140, 13143, 13144, 13148, 13149, 13150, 13154, 13155, 13159, 13160, 13161, 13162, 13166, 13167, 13171, 13172, 13173, 13174, 13178, 13179, 13183, 13184, 13185, 13186, 13189, 13190, 13191, 13192, 13198, 13199, 13200, 13201, 13205, 13206, 13207, 13208, 13209, 13210, 13218, 13219, 13221, 13222, 13228, 13230, 13231, 13232, 13234, 13235, 13239, 13240, 13241, 13242, 13246, 13247, 13251, 13252, 13253, 13258, 13259, 13260, 13261, 13264, 13265, 13266, 13270, 13271, 13272, 13273, 13274, 13275, 13276, 13279, 13281, 13283, 13285, 13292, 13293, 13296, 13298, 13299, 13300, 13303, 13304, 13306, 13311, 13313, 13314, 13315, 13316, 13320, 13321, 13322, 13323, 13324, 13328, 13329, 13332, 13334, 13337, 13340, 13341, 13342, 13344, 13348, 13353, 13354, 13355, 13356, 13357, 13362, 13363, 13364, 13370, 13371, 13372, 13377, 13378, 13379, 13382, 13384, 13386, 13387, 13388, 13390, 13391, 13397, 13401, 13402, 13403, 13404, 13405, 13410, 13412, 13416, 13417, 13418, 13419, 13422, 13424, 13425, 13427, 13431, 13433, 13435, 13438, 13439, 13440, 13441, 13442, 13443, 13444, 13445, 13448, 13454, 13458, 13460, 13464, 13466, 13467, 13469, 13471, 13472, 13473, 13477, 13479, 13480, 13482, 13483, 13484, 13485, 13489, 13490, 13495, 13496, 13498, 13499, 13500, 13506, 13507, 13508, 13509, 13516, 13517, 13518, 13519, 13520, 13522, 13523, 13528, 13529, 13530, 13531, 13533, 13540, 13541, 13542, 13545, 13546, 13547, 13552, 13553, 13554, 13557, 13558, 13559, 13564, 13565, 13569, 13570, 13571, 13572, 13574, 13575, 13579, 13580, 13581, 13582, 13583, 13588, 13589, 13593, 13594, 13595, 13596, 13600, 13601, 13605, 13606, 13607, 13608, 13610, 13616, 13617, 13618, 13620, 13621, 13625, 13626, 13627, 13630, 13631, 13633, 13637, 13639, 13643, 13645, 13646, 13649, 13650, 13651, 13654, 13655, 13656, 13657, 13659, 13661, 13662, 13663, 13664, 13666, 13669, 13670, 13672, 13673, 13674, 13675, 13678, 13679, 13680, 13685, 13686, 13687, 13690, 13691, 13692, 13697, 13698, 13699, 13701, 13702, 13708, 13709, 13710, 13711, 13712, 13715, 13716, 13717, 13718, 13719, 13724, 13725, 13726, 13732, 13734, 13735, 13736, 13737, 13742, 13743, 13746, 13747, 13749, 13750, 13751, 13753, 13754, 13755, 13758, 13760, 13761, 13762, 13763, 13764, 13766, 13769, 13771, 13772, 13773, 13775, 13778, 13781, 13782, 13783, 13784, 13788, 13791, 13792, 13793, 13797, 13798, 13802, 13803, 13804, 13805, 13808, 13809, 13810, 13813, 13814, 13815, 13817, 13818, 13819, 13820, 13825, 13827, 13828, 13830, 13834, 13835, 13836, 13840, 13843, 13844, 13845, 13846, 13848, 13850, 13853, 13855, 13857, 13859, 13860, 13861, 13865, 13870, 13871, 13876, 13877, 13878, 13879, 13880, 13884, 13885, 13886, 13891, 13892, 13893, 13894, 13896, 13900, 13902, 13903, 13906, 13908, 13911, 13913, 13918, 13919, 13920, 13922, 13925, 13926, 13933, 13935, 13936, 13939, 13940, 13941, 13944, 13946, 13948, 13949, 13950, 13952, 13954, 13959, 13961, 13963, 13965, 13966, 13967, 13968, 13970, 13971, 13973, 13974, 13975, 13977, 13982, 13985, 13986, 13987, 13989, 13990, 13994, 13995, 13997, 13998, 13999, 14000, 14002, 14003, 14008, 14011, 14012, 14013, 14017, 14018, 14022, 14023, 14024, 14025, 14026, 14030, 14031, 14035, 14036, 14037, 14038, 14040, 14047, 14048, 14052, 14053, 14054, 14055, 14062, 14064, 14066, 14068, 14070, 14071, 14072, 14073, 14075, 14076, 14078, 14079, 14081, 14083, 14084, 14086, 14089, 14093, 14094, 14095, 14097, 14102, 14103, 14104, 14108, 14110, 14112, 14113, 14115, 14116, 14120, 14121, 14124, 14125, 14129, 14132, 14133, 14134, 14136, 14138, 14141, 14143, 14144, 14145, 14147, 14148, 14150, 14154, 14155, 14156, 14159, 14160, 14161, 14163, 14167, 14168, 14169, 14173, 14176, 14177, 14178, 14179, 14183, 14186, 14187, 14188, 14189, 14190, 14191, 14194, 14200, 14201, 14202, 14203, 14204, 14205, 14207, 14211, 14212, 14213, 14215, 14217, 14218, 14219, 14220, 14224, 14225, 14226, 14228, 14231, 14232, 14233, 14234, 14235, 14239, 14241, 14243, 14247, 14248, 14254, 14259, 14260, 14261, 14262, 14263, 14269, 14271, 14272, 14273, 14277, 14278, 14283, 14284, 14287, 14288, 14291, 14296, 14298, 14300, 14301, 14302, 14303, 14305, 14313, 14314, 14315, 14317, 14320, 14321, 14322, 14323, 14324, 14325, 14326, 14331, 14333, 14334, 14336, 14337, 14340, 14341, 14342, 14346, 14347, 14349, 14350, 14353, 14355, 14356, 14357, 14360, 14361, 14362, 14364, 14369, 14370, 14371, 14375, 14376, 14377, 14384, 14385, 14386, 14387, 14388, 14390, 14394, 14400, 14401, 14405, 14406, 14407, 14408, 14410, 14411, 14413, 14415, 14418, 14422, 14424, 14427, 14428, 14429, 14430, 14431, 14438, 14440, 14441, 14442, 14444, 14446, 14449, 14451, 14452, 14456, 14460, 14465, 14467, 14468, 14470, 14476, 14477, 14478, 14481, 14483, 14485, 14488, 14489, 14491, 14493, 14495, 14498, 14501, 14502, 14504, 14506, 14508, 14509, 14511, 14512, 14513, 14514, 14516, 14521, 14523, 14525, 14527, 14529, 14531, 14532, 14533, 14535, 14536, 14537, 14538, 14541, 14544, 14546, 14547, 14548, 14549, 14550, 14551, 14552, 14557, 14558, 14559, 14560, 14566, 14567, 14568, 14569, 14570, 14571, 14576, 14577, 14581, 14584, 14586, 14589, 14590, 14591, 14593, 14595, 14596, 14597, 14599, 14602, 14603, 14606, 14607, 14609, 14611, 14612, 14614, 14615, 14616, 14617, 14620, 14622, 14623, 14625, 14630, 14631, 14633, 14634, 14635, 14636, 14638, 14639, 14641, 14643, 14645, 14646, 14648, 14651, 14653, 14654, 14655, 14656, 14661, 14663, 14664, 14666, 14669, 14670, 14671, 14674, 14675, 14679, 14682, 14684, 14685, 14687, 14691, 14693, 14695, 14700, 14701, 14702, 14703, 14705, 14706, 14707, 14710, 14712, 14714, 14715, 14716, 14717, 14720, 14721, 14722, 14724, 14726, 14729, 14730, 14734, 14737, 14738, 14739, 14740, 14741, 14742, 14745, 14746, 14749, 14752, 14754, 14756, 14757, 14759, 14762, 14764, 14766, 14768, 14769, 14773, 14774, 14776, 14779, 14782, 14783, 14784, 14785, 14786, 14788, 14789, 14793, 14795, 14796, 14799, 14801, 14802, 14804, 14808, 14809, 14811, 14814, 14816, 14818, 14819, 14820, 14821, 14824, 14825, 14827, 14828, 14830, 14831, 14835, 14836, 14837, 14838, 14839, 14842, 14844, 14845, 14847, 14849, 14850, 14852, 14853, 14856, 14860, 14861, 14862, 14864, 14866, 14869, 14872, 14874, 14876, 14878, 14880, 14882, 14884, 14885, 14886, 14888, 14889, 14891, 14893, 14899, 14901, 14902, 14905, 14906, 14908, 14909, 14911, 14912, 14915, 14917, 14920, 14924, 14928, 14932, 14934, 14936, 14939, 14943, 14946, 14949, 14951, 14955, 14959, 14963, 14965, 14967, 14972, 14976, 14980, 14983, 14987, 14991, 14997, 15000, 15004, 15008, 15010, 15014, 15016, 15020, 15025, 15026, 15027, 15028, 15029, 15030, 15031, 15033, 15035, 15037, 15039, 15040, 15042, 15046, 15047, 15048, 15050, 15051, 15052, 15054, 15055, 15056, 15057, 15060, 15061, 15065, 15067, 15069, 15071, 15073, 15075, 15077, 15079, 15080, 15081, 15083, 15086, 15089, 15091, 15092, 15094, 15096, 15097, 15099, 15101, 15105, 15106, 15110, 15112, 15113, 15115, 15116, 15118, 15119, 15123, 15125, 15126, 15127, 15129, 15131, 15133, 15136, 15137, 15139, 15142, 15144, 15146, 15149, 15151, 15152, 15155, 15158, 15159, 15163, 15166, 15167, 15169, 15171, 15174, 15176, 15178, 15180, 15183, 15185, 15186, 15188, 15190, 15195, 15199, 15202, 15204, 15207, 15210, 15212, 15214, 15216, 15217, 15219, 15222, 15225, 15227, 15229, 15232, 15233, 15234, 15236, 15239, 15240, 15241, 15242, 15246, 15248, 15249, 15254, 15255, 15258, 15261, 15262, 15263, 15265, 15266, 15268, 15272, 15276, 15278, 15279, 15280, 15281, 15282, 15284, 15286, 15287, 15288, 15289, 15290, 15291, 15293, 15295, 15298, 15299, 15300, 15301, 15302, 15303, 15307, 15308, 15312, 15313, 15314, 15317, 15318, 15319, 15324, 15325, 15326, 15328, 15332, 15333, 15334, 15338, 15339, 15340, 15341, 15345, 15346, 15350, 15351, 15352, 15353, 15357, 15358, 15362, 15363, 15364, 15365, 15369, 15372, 15373, 15374, 15376, 15377, 15381, 15382, 15383, 15384, 15387, 15391, 15392, 15394, 15395, 15397, 15402, 15404, 15406, 15407, 15408, 15409, 15410, 15412, 15413, 15414, 15417, 15418, 15419, 15420, 15423, 15424, 15426, 15428, 15429, 15432, 15434, 15435, 15436, 15438, 15445, 15446, 15448, 15453, 15454, 15455, 15457, 15458, 15461, 15464, 15465, 15468, 15469, 15470, 15471, 15472, 15473, 15475, 15476, 15479, 15482, 15484, 15485, 15487, 15488, 15489, 15490, 15493, 15494, 15495, 15496, 15497, 15498, 15501, 15503, 15506, 15508, 15510, 15514, 15516, 15517, 15519, 15520, 15521, 15523, 15524, 15525, 15528, 15532, 15538, 15540, 15541, 15542, 15543, 15546, 15547, 15548, 15549, 15553, 15556, 15557, 15558, 15559, 15560, 15561, 15563, 15564, 15565, 15567, 15568, 15572, 15573, 15574, 15575, 15576, 15580, 15583, 15584, 15585, 15586, 15594, 15596, 15598, 15600, 15601, 15602, 15603, 15606, 15607, 15609, 15613, 15614, 15615, 15617, 15622, 15625, 15628, 15629, 15630, 15634, 15635, 15640, 15641, 15643, 15644, 15647, 15648, 15649, 15652, 15653, 15655, 15657, 15661, 15663, 15664, 15665, 15672, 15674, 15675, 15676, 15677, 15681, 15682, 15683, 15684, 15688, 15692, 15694, 15695, 15696, 15697, 15699, 15701, 15702, 15704, 15706, 15710, 15711, 15715, 15717, 15719, 15723, 15724, 15725, 15726, 15730, 15732, 15733, 15734, 15735, 15737, 15738, 15742, 15743, 15744, 15745, 15746, 15751, 15752, 15753, 15755, 15756, 15763, 15764, 15765, 15768, 15771, 15772, 15773, 15774, 15775, 15777, 15778, 15779, 15781, 15782, 15783, 15785, 15786, 15787, 15789, 15791, 15793, 15795, 15796, 15799, 15802, 15803, 15804, 15808, 15809, 15810, 15813, 15814, 15815, 15816, 15817, 15818, 15822, 15823, 15825, 15826, 15827, 15831, 15833, 15835, 15836, 15837, 15840, 15846, 15847, 15848, 15849, 15851, 15856, 15857, 15858, 15859, 15861, 15863, 15869, 15870, 15872, 15873, 15876, 15878, 15879, 15885, 15888, 15889, 15892, 15893, 15896, 15897, 15899, 15916, 15920, 15922, 15924, 15925, 15928, 15930, 15931, 15935, 15936, 15938, 15941, 15943, 15944, 15946, 15953, 15954, 15955, 15957, 15958, 15959, 15961, 15963, 15965, 15967, 15968, 15970, 15971, 15975, 15977, 15979, 15982, 15984, 15986, 15992, 15996, 16000, 16004, 16008, 16012, 16016, 16022, 16025, 16026, 16027, 16028, 16029, 16030, 16032, 16033, 16034, 16036, 16037, 16040, 16044, 16047, 16051, 16053, 16056, 16060, 16062, 16066, 16070, 16074, 16075, 16076, 16079, 16080, 16081, 16083, 16084, 16085, 16087, 16089, 16090, 16091, 16093, 16094, 16098, 16099, 16101, 16103, 16104, 16105, 16107, 16108, 16109, 16111, 16112, 16113, 16114, 16116, 16118, 16121, 16122, 16124, 16125, 16127, 16129, 16130, 16131, 16134, 16138, 16139, 16141, 16143, 16146, 16148, 16151, 16152, 16153, 16155, 16156, 16157, 16160, 16161, 16162, 16164, 16165, 16168, 16170, 16171, 16172, 16174, 16175, 16176, 16179, 16181, 16184, 16185, 16187, 16188, 16189, 16192, 16193, 16195, 16198, 16200, 16201, 16205, 16209, 16210, 16212, 16213, 16215, 16216, 16217, 16219, 16220, 16223, 16224, 16226, 16228, 16229, 16230, 16232, 16233, 16234, 16236, 16238, 16240, 16241, 16242, 16244, 16245, 16246, 16248, 16249, 16251, 16252, 16254, 16255, 16257, 16261, 16262, 16265, 16267, 16269, 16270, 16272, 16273, 16275, 16278, 16279, 16280, 16282, 16283, 16284, 16286, 16289, 16292, 16294, 16296, 16297, 16300, 16301, 16304, 16306, 16308, 16310, 16313, 16315, 16316, 16319, 16322, 16326, 16327, 16329, 16332, 16337, 16338, 16341, 16345, 16346, 16348, 16350, 16354, 16356, 16359, 16360, 16361, 16363, 16365, 16368, 16370, 16372, 16375, 16376, 16377, 16378, 16380, 16382, 16383, 16387, 16389, 16390, 16392, 16393, 16395, 16397, 16399, 16401, 16405, 16408, 16409, 16412, 16414, 16416, 16418, 16419, 16422, 16424, 16425, 16426, 16429, 16431, 16432, 16434, 16435, 16438, 16444, 16445, 16446, 16450, 16451, 16452, 16456, 16457, 16461, 16462, 16463, 16464, 16468, 16469, 16470, 16471, 16472, 16473, 16476, 16477, 16478, 16483, 16484, 16485, 16487, 16488, 16489, 16491, 16492, 16494, 16495, 16496, 16497, 16498, 16501, 16502, 16503, 16507, 16508, 16509, 16510, 16513, 16517, 16520, 16521, 16524, 16525, 16526, 16527, 16531, 16532, 16534, 16535, 16536, 16540, 16541, 16542, 16544, 16548, 16549, 16550, 16551, 16555, 16556, 16561, 16562, 16563, 16564, 16566, 16573, 16574, 16575, 16576, 16577, 16579, 16580, 16582, 16584, 16588, 16591, 16592, 16593, 16595, 16596, 16603, 16605, 16612, 16613, 16614, 16615, 16616, 16620, 16621, 16622, 16623, 16624, 16625, 16630, 16631, 16634, 16636, 16637, 16638, 16641, 16646, 16648, 16649, 16650, 16651, 16653, 16654, 16655, 16656, 16657, 16659, 16664, 16665, 16667, 16668, 16672, 16673, 16674, 16679, 16680, 16681, 16683, 16684, 16686, 16687, 16690, 16691, 16693, 16696, 16699, 16700, 16701, 16706, 16707, 16708, 16709, 16711, 16718, 16719, 16721, 16722, 16723, 16729, 16730, 16731, 16732, 16736, 16737, 16739, 16740, 16745, 16746, 16747, 16748, 16752, 16754, 16755, 16757, 16760, 16762, 16767, 16769, 16771, 16773, 16774, 16777, 16778, 16781, 16782, 16785, 16786, 16788, 16789, 16791, 16792, 16793, 16795, 16796, 16797, 16799, 16804, 16807, 16808, 16809, 16810, 16812, 16814, 16815, 16816, 16818, 16820, 16826, 16827, 16832, 16833, 16838, 16839, 16841, 16842, 16847, 16848, 16849, 16853, 16856, 16857, 16858, 16859, 16860, 16861, 16865, 16866, 16869, 16871, 16872, 16873, 16874, 16880, 16881, 16884, 16889, 16890, 16892, 16894, 16898, 16899, 16900, 16901, 16902, 16903, 16905, 16909, 16910, 16911, 16912, 16914, 16917, 16918, 16920, 16921, 16922, 16925, 16927, 16930, 16933, 16934, 16935, 16938, 16940, 16941, 16942, 16944, 16947, 16948, 16951, 16953, 16954, 16960, 16961, 16962, 16963, 16965, 16966, 16968, 16970, 16974, 16977, 16979, 16982, 16983, 16987, 16989, 16991, 16993, 16994, 16995, 17000, 17001, 17004, 17005, 17007, 17011, 17012, 17013, 17017, 17018, 17022, 17023, 17026, 17030, 17032, 17033, 17035, 17036, 17038, 17041, 17046, 17047, 17051, 17055, 17059, 17060, 17062, 17064, 17066, 17067, 17069, 17076, 17077, 17078, 17082, 17083, 17084, 17086, 17089, 17090, 17092, 17093, 17096, 17098, 17099, 17101, 17103, 17104, 17106, 17107, 17109, 17110, 17111, 17114, 17117, 17123, 17125, 17126, 17128, 17131, 17132, 17138, 17141, 17143, 17144, 17145, 17147, 17149, 17151, 17152, 17154, 17155, 17157, 17159, 17162, 17163, 17165, 17167, 17168, 17169, 17171, 17173, 17175, 17177, 17179, 17181, 17184, 17186, 17189, 17191, 17194, 17196, 17199, 17200, 17203, 17204, 17206, 17207, 17209, 17210, 17211, 17214, 17216, 17219, 17220, 17222, 17225, 17227, 17228, 17229, 17230, 17232, 17233, 17234, 17236, 17237, 17238, 17239, 17240, 17242, 17243, 17244, 17248, 17252, 17254, 17256, 17258, 17259, 17262, 17267, 17271, 17275, 17277, 17279, 17280, 17281, 17282, 17283, 17285, 17286, 17287, 17289, 17290, 17291, 17293, 17295, 17296, 17297, 17299, 17300, 17301, 17303, 17304, 17307, 17313, 17317, 17321, 17323, 17327, 17329, 17330, 17336, 17341, 17342, 17343, 17345, 17346, 17348, 17350, 17352, 17353, 17355, 17357, 17358, 17363, 17364, 17366, 17368, 17370, 17376, 17380, 17382, 17384, 17386, 17387, 17389, 17391, 17393, 17394, 17395, 17397, 17400, 17403, 17405, 17406, 17407, 17409, 17412, 17415, 17418, 17420, 17421, 17424, 17426, 17428, 17432, 17434, 17435, 17436, 17438, 17440, 17442, 17443, 17444, 17446, 17447, 17449, 17451, 17453, 17454, 17458, 17459, 17460, 17461, 17462, 17463, 17466, 17468, 17470, 17472, 17473, 17475, 17476, 17478, 17480, 17482, 17484, 17488, 17489, 17491, 17493, 17494, 17495, 17496, 17498, 17500, 17502, 17503, 17504, 17507, 17508, 17509, 17510, 17512, 17515, 17517, 17519, 17520, 17521, 17523, 17525, 17528, 17530, 17531, 17533, 17536, 17539, 17541, 17543, 17545, 17547, 17549, 17552, 17554, 17555, 17558, 17560, 17562, 17563, 17566, 17569, 17572, 17575, 17578, 17579, 17582, 17585, 17588, 17590, 17592, 17594, 17596, 17597, 17600, 17603, 17610, 17611, 17612, 17613, 17615, 17617, 17619, 17621, 17623, 17625, 17626, 17627, 17629, 17631, 17632, 17633, 17636, 17638, 17641, 17644, 17647, 17648, 17650, 17654, 17656, 17657, 17659, 17661, 17663, 17664, 17666, 17668, 17669, 17671, 17672, 17674, 17675, 17678, 17680, 17681, 17683, 17686, 17688, 17689, 17691, 17693, 17695, 17696, 17699, 17700, 17701, 17702, 17704, 17705, 17706, 17708, 17710, 17713, 17714, 17716, 17718, 17720, 17723, 17725, 17726, 17728, 17731, 17733, 17734, 17736, 17737, 17739, 17742, 17744, 17746, 17748, 17749, 17750, 17752, 17754, 17756, 17759, 17761, 17763, 17764, 17765, 17766, 17768, 17769, 17770, 17772, 17775, 17777, 17778, 17780, 17783, 17785, 17786, 17788, 17789, 17790, 17791, 17795, 17799, 17801, 17803, 17805, 17806, 17808, 17810, 17811, 17813, 17816, 17818, 17819, 17821, 17823, 17824, 17825, 17827, 17830, 17832, 17833, 17835, 17836, 17837, 17840, 17841, 17843, 17845, 17846, 17848, 17850, 17851, 17852, 17854, 17856, 17858, 17859, 17863, 17867, 17871, 17873, 17877, 17878, 17881, 17884, 17886, 17888, 17890, 17893, 17895, 17897, 17899, 17902, 17905, 17906, 17909, 17910, 17912, 17914, 17915, 17917, 17919, 17920, 17923, 17924, 17926, 17928, 17931, 17932, 17935, 17938, 17941, 17943, 17944, 17948, 17950, 17952, 17954, 17955, 17957, 17958, 17959, 17962, 17963, 17964, 17966, 17969, 17971, 17973, 17975, 17977, 17978, 17980, 17983, 17985, 17986, 17988, 17991, 17993, 17994, 17996, 17997, 17998, 18000, 18004, 18006, 18007, 18009, 18010, 18012, 18013, 18016, 18017, 18019, 18020, 18022, 18025, 18027, 18028, 18030, 18032, 18034, 18039, 18040, 18041, 18042, 18043, 18044, 18045, 18048, 18050, 18052, 18053, 18054, 18056, 18062, 18065, 18066, 18068, 18069, 18071, 18074, 18078, 18079, 18083, 18085, 18086, 18087, 18088, 18089, 18094, 18095, 18096, 18097, 18100, 18102, 18109, 18112, 18113, 18118, 18119, 18120, 18122, 18123, 18124, 18128, 18129, 18130, 18134, 18135, 18136, 18139, 18141, 18146, 18147, 18149, 18150, 18151, 18152, 18154, 18155, 18160, 18164, 18165, 18167, 18169, 18173, 18174, 18175, 18179, 18181, 18182, 18183, 18184, 18189, 18196, 18197, 18198, 18201, 18202, 18203, 18204, 18208, 18209, 18213, 18214, 18215, 18216, 18220, 18221, 18225, 18226, 18227, 18228, 18231, 18232, 18233, 18234, 18238, 18239, 18243, 18244, 18245, 18250, 18251, 18252, 18256, 18257, 18261, 18262, 18263, 18264, 18269, 18270, 18271, 18275, 18278, 18279, 18280, 18284, 18285, 18286, 18290, 18292, 18293, 18296, 18297, 18298, 18299, 18306, 18307, 18309, 18312, 18313, 18314, 18324, 18325, 18326, 18327, 18330, 18331, 18333, 18334, 18335, 18338, 18340, 18341, 18345, 18346, 18347, 18350, 18351, 18352, 18357, 18358, 18359, 18360, 18361, 18365, 18369, 18370, 18371, 18372, 18373, 18374, 18377, 18379, 18380, 18385, 18386, 18388, 18389, 18394, 18396, 18397, 18398, 18399, 18400, 18402, 18403, 18411, 18412, 18413, 18415, 18416, 18417, 18422, 18423, 18427, 18429, 18430, 18431, 18433, 18434, 18441, 18442, 18443, 18444, 18445, 18447, 18449, 18451, 18452, 18453, 18457, 18458, 18459, 18460, 18461, 18467, 18472, 18473, 18474, 18479, 18480, 18481, 18482, 18484, 18489, 18490, 18491, 18492, 18493, 18496, 18497, 18498, 18500, 18503, 18505, 18513, 18514, 18515, 18518, 18519, 18520, 18526, 18527, 18531, 18532, 18534, 18535, 18540, 18541, 18542, 18543, 18544, 18547, 18549, 18550, 18551, 18552, 18555, 18556, 18557, 18560, 18562, 18567, 18568, 18569, 18572, 18573, 18574, 18585, 18586, 18587, 18591, 18592, 18593, 18594, 18595, 18599, 18602, 18603, 18607, 18608, 18609, 18615, 18616, 18617, 18618, 18620, 18621, 18622, 18623, 18626, 18627, 18628, 18629, 18633, 18636, 18637, 18638, 18639, 18640, 18641, 18645, 18646, 18650, 18651, 18652, 18653, 18657, 18658, 18662, 18663, 18664, 18669, 18670, 18671, 18672, 18673, 18677, 18678, 18679, 18684, 18685, 18686, 18687, 18689, 18690, 18691, 18699, 18700, 18701, 18703, 18706, 18708, 18709, 18710, 18712, 18713, 18718, 18722, 18725, 18726, 18727, 18728, 18732, 18733, 18736, 18738, 18739, 18741, 18742, 18748, 18750, 18751, 18752, 18754, 18756, 18758, 18759, 18763, 18767, 18770, 18771, 18772, 18774, 18777, 18778, 18783, 18785, 18788, 18790, 18793, 18795, 18801, 18803, 18804, 18810, 18811, 18812, 18814, 18815, 18817, 18818, 18820, 18821, 18824, 18826, 18827, 18829, 18834, 18835, 18836, 18838, 18841, 18843, 18845, 18846, 18849, 18850, 18851, 18852, 18857, 18859, 18861, 18862, 18863, 18864, 18865, 18867, 18872, 18873, 18876, 18878, 18880, 18886, 18887, 18891, 18892, 18893, 18895, 18897, 18903, 18904, 18905, 18906, 18911, 18912, 18913, 18916, 18917, 18918, 18923, 18924, 18925, 18928, 18929, 18930, 18931, 18935, 18936, 18940, 18941, 18942, 18943, 18947, 18948, 18952, 18953, 18954, 18955, 18959, 18960, 18965, 18966, 18967, 18970, 18971, 18972, 18973, 18977, 18978, 18982, 18983, 18984, 18985, 18989, 18990, 18991, 18996, 18997, 19001, 19002, 19003, 19004, 19008, 19010, 19011, 19012, 19018, 19020, 19021, 19023, 19026, 19027, 19031, 19035, 19037, 19038, 19039, 19040, 19041, 19048, 19050, 19051, 19053, 19054, 19055, 19060, 19063, 19064, 19065, 19067, 19068, 19070, 19072, 19073, 19077, 19078, 19079, 19080, 19086, 19087, 19088, 19089, 19093, 19094, 19098, 19099, 19100, 19101, 19104, 19106, 19108, 19109, 19110, 19112, 19114, 19116, 19117, 19118, 19119, 19125, 19126, 19127, 19128, 19132, 19133, 19134, 19141, 19142, 19143, 19144, 19145, 19147, 19151, 19153, 19154, 19155, 19156, 19157, 19158, 19162, 19167, 19172, 19173, 19176, 19178, 19179, 19180, 19184, 19187, 19188, 19189, 19190, 19191, 19192, 19196, 19197, 19202, 19203, 19204, 19207, 19211, 19212, 19214, 19215, 19216, 19218, 19219, 19222, 19224, 19225, 19230, 19232, 19236, 19237, 19238, 19241, 19242, 19243, 19248, 19249, 19250, 19254, 19255, 19256, 19257, 19261, 19262, 19267, 19268, 19270, 19273, 19275, 19276, 19277, 19278, 19282, 19283, 19286, 19288, 19289, 19290, 19294, 19295, 19296, 19297, 19301, 19302, 19303, 19308, 19309, 19310, 19313, 19314, 19315, 19320, 19321, 19323, 19325, 19330, 19331, 19332, 19333, 19338, 19339, 19341, 19345, 19346, 19347, 19348, 19352, 19353, 19354, 19360, 19363, 19364, 19366, 19367, 19368, 19369, 19373, 19374, 19378, 19379, 19380, 19381, 19385, 19386, 19391, 19392, 19393, 19396, 19397, 19398, 19403, 19404, 19405, 19406, 19407, 19408, 19410, 19411, 19415, 19416, 19417, 19418, 19419, 19420, 19421, 19422, 19424, 19425, 19428, 19430, 19431, 19432, 19436, 19437, 19438, 19442, 19445, 19446, 19447, 19456, 19460, 19462, 19463, 19464, 19465, 19466, 19467, 19469, 19470, 19473, 19475, 19476, 19478, 19479, 19482, 19483, 19484, 19485, 19488, 19489, 19491, 19493, 19496, 19499, 19500, 19501, 19505, 19507, 19509, 19513, 19514, 19515, 19516, 19517, 19518, 19522, 19524, 19525, 19529, 19530, 19533, 19536, 19538, 19541, 19544, 19545, 19546, 19548, 19552, 19553, 19555, 19556, 19557, 19559, 19560, 19561, 19562, 19565, 19567, 19572, 19577, 19579, 19580, 19581, 19583, 19584, 19585, 19586, 19589, 19592, 19594, 19599, 19603, 19604, 19612, 19613, 19614, 19615, 19617, 19618, 19619, 19620, 19621, 19622, 19628, 19630, 19633, 19636, 19637, 19638, 19641, 19643, 19644, 19645, 19648, 19649, 19651, 19654, 19655, 19656, 19657, 19658, 19659, 19664, 19666, 19669, 19670, 19674, 19675, 19676, 19677, 19679, 19680, 19683, 19684, 19686, 19689, 19690, 19692, 19693, 19695, 19696, 19697, 19698, 19702, 19703, 19704, 19707, 19708, 19709, 19714, 19715, 19716, 19718, 19719, 19723, 19725, 19726, 19727, 19728, 19734, 19735, 19736, 19737, 19741, 19742, 19743, 19746, 19747, 19748, 19753, 19754, 19755, 19756, 19757, 19758, 19760, 19761, 19765, 19766, 19771, 19772, 19773, 19776, 19777, 19779, 19781, 19783, 19784, 19785, 19790, 19791, 19798, 19799, 19800, 19801, 19802, 19803, 19804, 19806, 19812, 19819, 19820, 19821, 19822, 19824, 19829, 19830, 19832, 19833, 19834, 19835, 19836, 19837, 19838, 19843, 19844, 19845, 19847, 19852, 19853, 19854, 19857, 19858, 19861, 19862, 19865, 19866, 19868, 19869, 19871, 19872, 19874, 19879, 19880, 19881, 19885, 19888, 19889, 19890, 19892, 19895, 19897, 19899, 19902, 19904, 19906, 19910, 19911, 19914, 19916, 19917, 19918, 19919, 19920, 19922, 19924, 19925, 19927, 19932, 19933, 19937, 19941, 19942, 19943, 19944, 19945, 19946, 19947, 19949, 19954, 19955, 19956, 19959, 19964, 19965, 19967, 19968, 19970, 19972, 19975, 19976, 19977, 19978, 19979, 19980, 19981, 19983, 19986, 19988, 19989, 19990, 19993, 19995, 19996, 19997, 19998, 20000, 20002, 20004, 20007, 20008, 20009, 20011, 20012, 20013, 20015, 20017, 20021, 20022, 20024, 20027, 20029, 20032, 20033, 20035, 20037, 20039, 20040, 20042, 20043, 20044, 20046, 20049, 20050, 20051, 20052, 20053, 20055, 20058, 20059, 20060, 20062, 20066, 20068, 20071, 20073, 20074, 20076, 20078, 20079, 20080, 20081, 20082, 20083, 20087, 20089, 20090, 20091, 20093, 20095, 20097, 20106, 20111, 20113, 20117, 20120, 20121, 20123, 20124, 20125, 20126, 20128, 20129, 20130, 20131, 20132, 20133, 20136, 20137, 20138, 20139, 20140, 20142, 20144, 20146, 20148, 20150, 20153, 20155, 20157, 20159, 20160, 20162, 20163, 20165, 20166, 20168, 20169, 20170, 20174, 20177, 20178, 20179, 20180, 20183, 20189, 20190, 20194, 20195, 20197, 20198, 20199, 20201, 20202, 20203, 20204, 20208, 20209, 20210, 20211, 20215, 20216, 20218, 20219, 20221, 20223, 20224, 20225, 20228, 20229, 20230, 20235, 20237, 20240, 20241, 20243, 20245, 20246, 20247, 20248, 20252, 20254, 20255, 20256, 20259, 20260, 20262, 20263, 20264, 20266, 20269, 20270, 20272, 20273, 20274, 20276, 20279, 20281, 20286, 20287, 20289, 20292, 20297, 20298, 20301, 20302, 20303, 20306, 20307, 20308, 20310, 20311, 20315, 20316, 20318, 20319, 20321, 20322, 20325, 20327, 20329, 20331, 20332, 20335, 20336, 20337, 20338, 20339, 20340, 20341, 20342, 20343, 20344, 20345, 20346, 20349, 20351, 20352, 20353, 20354, 20355, 20357, 20358, 20359, 20360, 20361, 20362, 20363, 20364, 20365, 20366, 20367, 20369, 20370, 20371, 20372, 20373, 20374, 20375, 20376, 20377, 20378, 20379, 20380, 20381, 20382, 20383, 20384, 20385, 20386, 20387, 20388, 20389, 20390, 20391, 20392, 20393, 20394, 20395, 20396, 20397, 20398, 20399, 20400, 20401, 20402, 20403, 20404, 20405, 20406, 20407, 20408, 20409, 20410, 20411, 20412, 20413, 20414, 20415, 20416, 20417, 20418, 20419, 20420, 20421, 20422, 20423, 20424, 20425, 20426, 20456, 20457, 20458, 20459, 20460, 20461, 20462, 20463, 20464, 20465, 20466, 20467, 20468, 20475, 20476, 20477, 20478, 20479, 20480, 20481, 20482, 20483, 20484, 20485, 20486, 20487, 20488, 20489, 20490, 20491, 20492, 20503, 20504, 20505, 20506, 20507, 20508, 20509, 20510, 20511, 20512, 20513, 20514, 20515, 20516, 20517, 20518, 20519, 20520, 20521, 20522, 20523, 20524, 20525, 20526, 20527, 20528, 20529, 20530, 20531, 20532, 20533, 20534, 20536, 20537, 20538, 20539, 20540, 20541, 20542, 20543, 20544, 20545, 20546, 20547, 20548, 20549, 20550, 20551, 20552, 20553, 20554, 20555, 20556, 20557, 20558, 20559, 20560, 20561, 20562, 20563, 20564, 20565, 20566, 20567, 20568, 20569, 20570, 20571, 20572, 20574, 20575, 20576, 20577, 20578, 20579, 20580, 20581, 20582, 20583, 20584, 20585, 20586, 20587, 20588, 20589, 20590, 20591, 20592, 20593, 20594, 20595, 20596, 20597, 20598, 20599, 20600, 20601, 20602, 20603, 20604, 20605, 20606, 20607, 20608, 20609, 20610, 20611, 20613, 20616, 20619, 20622, 20625, 20628, 20631, 20634, 20637, 20640, 20656, 20657, 20658, 20659, 20660, 20661, 

info: 1 property vectors copied, 0 vectors skipped.
info: 17 hole(s) detected within the mesh
info: Remark: Nodes n0 and n1 are considered to be 'collapsible' if the dist(n0, n1) is less than minimal edge length of elements in the mesh times 1e-6

Boundary and materials extraction using OpenGeoSys tools

ExtractBoundary Extracts all outer surface faces from a 3D mesh and creates a boundary-only mesh. Useful for defining where boundary conditions will be applied.
RemoveMeshElements Removes elements based on coordinate filters (e.g., x-min, x-max). Used to isolate specific boundary patches for inflow/outflow or other localized conditions.
ExtractMaterials Extracts submeshes from a mesh by selecting elements with specific material IDs, enabling features such as boreholes, fractures, or layers to be used for boundary conditions.

domain_size = 100.0… (click to toggle)

domain_size = 100.0
mins = np.array([0.0, 0.0, 0.0])
maxs = np.array([domain_size, domain_size, domain_size])
n_fractures = 10


boundary_tol_fraction = 1e-6
boundary_tol = boundary_tol_fraction * domain_size
x_outlet = maxs[0] - boundary_tol
x_inlet = mins[0] + boundary_tol

ot.cli().ExtractBoundary(i="mixed_dimensional_grid_2.vtu", o="boundaries.vtu")… (click to toggle)

ot.cli().ExtractBoundary(i="mixed_dimensional_grid_2.vtu", o="boundaries.vtu")
ot.cli().removeMeshElements(i="boundaries.vtu", o="x_inlet.vtu", **{"x-min": x_inlet})
ot.cli().removeMeshElements(
    i="boundaries.vtu", o="x_outlet.vtu", **{"x-min": x_outlet, "invert": True}
)

info: Mesh read: 10340 nodes, 20662 elements.
info: 1 property vectors copied, 0 vectors skipped.
info: Created surface mesh: 1132 nodes, 1154 elements.
info: Mesh read: 1132 nodes, 1154 elements.
info: Bounding box of "boundaries" is
x = [-0.000000,100.000000]
y = [0.000000,100.000000]
z = [0.000000,100.000000]
info: 988 elements found.
info: Removing total 988 elements...
info: 166 elements remain in mesh.
info: Removing total 967 nodes...
info: Mesh read: 1132 nodes, 1154 elements.
info: Bounding box of "boundaries" is
x = [-0.000000,100.000000]
y = [0.000000,100.000000]
z = [0.000000,100.000000]
info: 1129 elements found.
info: Removing total 1129 elements...
info: 25 elements remain in mesh.
info: Removing total 1105 nodes...

ot.cli().ExtractMaterials(… (click to toggle)

ot.cli().ExtractMaterials(
    i="mixed_dimensional_grid_2.vtu", o="borehole_right.vtu", m=n_fractures
)
ot.cli().ExtractMaterials(
    i="mixed_dimensional_grid_2.vtu", o="borehole_left.vtu", m=n_fractures + 1
)
ot.cli().identifySubdomains(
    "-m",
    "mixed_dimensional_grid_2.vtu",
    "-s",
    str(1e-10),
    "--",
    f"borehole_left_Layer{n_fractures+1}.vtu",
    f"borehole_right_Layer{n_fractures}.vtu",
)

info: Extracting material group 10...
info: Removing total 20490 elements...
info: 172 elements remain in mesh.
info: Removing total 10192 nodes...
info: Extracting material group 11...
info: Removing total 20512 elements...
info: 150 elements remain in mesh.
info: Removing total 10206 nodes...

info: Mesh reading time: 0.0164842 s
info: MeshNodeSearcher construction time: 0.000560847 s
info: identifySubdomainMesh(): identifySubdomainMeshNodes took 2.8122e-05 s
info: identifySubdomainMesh(): identifySubdomainMeshElements took 0.000886289 s
info: identifySubdomainMesh(): identifySubdomainMeshNodes took 2.2213e-05 s
info: identifySubdomainMesh(): identifySubdomainMeshElements took 0.000462168 s
info: identifySubdomains time: 0.00151786 s
info: writing time: 0.00186096 s
info: Entire run time: 0.0224692 s

Plot DFN and boundaries

os.chdir(orig_dir)… (click to toggle)

os.chdir(orig_dir)
try:
    pv.set_jupyter_backend("static")
except Exception as e:
    print("PyVista backend not set:", e)

DFN_2D = pv.read(f"{out_dir}/mixed_dimensional_grid_2.vtu")
x_inlet = pv.read(f"{out_dir}/x_inlet.vtu")
x_outlet = pv.read(f"{out_dir}/x_outlet.vtu")

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(
    DFN_2D,
    scalars="MaterialIDs",
    cmap="tab20",
    opacity=0.4,
    scalar_bar_args={"title": "MaterialIDs", "vertical": True},
)
plotter.add_mesh(DFN_2D.outline(), color="black", line_width=2)
plotter.add_mesh(
    x_inlet,
    color="blue",
    opacity=1.0,
    line_width=5,
    render_lines_as_tubes=True,
)
plotter.add_mesh(
    x_outlet,
    color="red",
    opacity=1.0,
    line_width=5,
    render_lines_as_tubes=True,
)

plotter.show_axes()
plotter.enable_parallel_projection()
plotter.view_isometric()
output_path = out_dir.joinpath("Concentration.png")
plotter.screenshot(str(output_path))
plotter.show()
plotter.close()

[0m[33m2026-05-18 15:52:19.897 (   3.422s) [    7FA918613400]vtkXOpenGLRenderWindow.:1460  WARN| bad X server connection. DISPLAY=[0m

png

Generate and assign random fracture properties with Gaussian distributions

We assign random fracture properties (e.g., width and permeability) by calculating each mean value, $\mu$, from one of three user-defined distributions: Uniform, clipped Gaussian, or truncated Gaussian, constrained within specified bounds. The standard deviation is then calculated as $\sigma = \mathrm{rel_{std}} \times \mu$, with $\mathrm{rel_{std}}$ fixed at 30%, and cell values are subsequently sampled from a Gaussian distribution $\bigl(\mu, \sigma\bigr)$ and mapped to the mesh. Permeability can be configured as either isotropic or anisotropic.

Remark: Random fields with spatial correlation would be a more suitable representation of fracture heterogeneity. In the context of OGS see [8,9].

def sample_truncnorm(a, b, mu, sigma, rng, size=None):… (click to toggle)

def sample_truncnorm(a, b, mu, sigma, rng, size=None):
    if size is None:
        # scalar
        while True:
            x = rng.normal(mu, sigma)
            if a <= x <= b:
                return x
    size = int(size)
    out = np.empty(size)
    filled = 0
    batch = max(1024, size)
    while filled < size:
        x = rng.normal(mu, sigma, batch)
        x = x[(x >= a) & (x <= b)]
        take = min(size - filled, x.size)
        if take:
            out[filled : filled + take] = x[:take]
            filled += take
    return out


def sample_mean(a, b, rng, method="uniform"):
    if method == "uniform":
        return rng.uniform(a, b)
    if method == "constant":
        return 0.5 * (a + b)

    mu = 0.5 * (a + b)
    sigma = (b - a) / 6.0

    if method == "gaussian":
        return np.clip(rng.normal(mu, sigma), a, b)
    if method == "truncated":
        return sample_truncnorm(a, b, mu, sigma, rng)

    msg = "Unknown sampling method"
    raise ValueError(msg)


def generate_random_fracture_stats(
    fracture_ids,
    width_range=(5e-6, 1e-5),
    k_range=(1e-14, 1e-12),
    rel_std=0.3,
    mean_dist="uniform",
    rng=None,
):
    rng = rng or default_rng()
    stats = {}
    for mat_id in fracture_ids:
        w_min, w_max = width_range
        width_mean = sample_mean(w_min, w_max, rng, method=mean_dist)
        k_min, k_max = k_range
        k_mean = sample_mean(k_min, k_max, rng, method=mean_dist)
        stats[mat_id] = {
            "width_mean": width_mean,
            "width_std": rel_std * width_mean,
            "width_range": (w_min, w_max),
            "k_mean": k_mean,
            "k_std": rel_std * k_mean,
            "k_range": (k_min, k_max),
            "method": mean_dist,
        }
    return stats


def assign_cellwise_random_props(mesh, stats_dict, clip_min=1e-20, seed=42):
    rng = default_rng(seed)
    MaterialIDs = mesh.cell_data["MaterialIDs"]
    n_cells = len(MaterialIDs)
    width_ic = np.zeros(n_cells)
    permeability = np.zeros(n_cells)

    for mat_id, stats in stats_dict.items():
        idx = np.where(MaterialIDs == mat_id)[0]
        n = len(idx)
        method = stats.get("method", "uniform")
        w_min, w_max = stats["width_range"]
        w_mean, w_std = stats["width_mean"], stats["width_std"]
        k_min, k_max = stats["k_range"]
        k_mean, k_std = stats["k_mean"], stats["k_std"]

        if method == "uniform":
            width_vals = rng.uniform(w_min, w_max, n)
            k_vals = rng.uniform(k_min, k_max, n)
        elif method == "constant":
            width_vals = np.full(n, w_mean)
            k_vals = np.full(n, k_mean)
        elif method == "truncated":
            width_vals = sample_truncnorm(w_min, w_max, w_mean, w_std, rng, size=n)
            k_vals = sample_truncnorm(k_min, k_max, k_mean, k_std, rng, size=n)
        elif method == "gaussian":
            width_vals = np.clip(rng.normal(w_mean, w_std, n), w_min, w_max)
            k_vals = np.clip(rng.normal(k_mean, k_std, n), k_min, k_max)
        else:
            msg = "Unknown sampling method"
            raise ValueError(msg)

        width_ic[idx] = np.clip(width_vals, clip_min, None)
        permeability[idx] = k_vals

    mesh.cell_data["width_ic"] = width_ic
    mesh.cell_data["permeability_ic"] = permeability

input_mesh = "mixed_dimensional_grid_2.vtu"… (click to toggle)

input_mesh = "mixed_dimensional_grid_2.vtu"
output_mesh = "mixed_dimensional_grid_2_updated.vtu"

mesh = pv.read(Path(out_dir, input_mesh))

material_ids = np.unique(mesh.cell_data["MaterialIDs"])
fracture_ids = list(material_ids)

SEED = 20250818… (click to toggle)

SEED = 20250818
rng = default_rng(SEED)

fracture_stats = generate_random_fracture_stats(
    fracture_ids,
    width_range=(1e-6, 1e-5),
    k_range=(1e-16, 1e-14),
    rel_std=0.3,  # standard deviation is 30 % of mean value.
    mean_dist="gaussian",  # constant, uniform, gaussian, or truncated
    rng=rng,
)

assign_cellwise_random_props(mesh, fracture_stats, seed=SEED)
mesh.save(Path(out_dir, output_mesh), binary=False)

base_scalar_bar_args = {… (click to toggle)

base_scalar_bar_args = {
    "vertical": True,
    "position_x": 0.9,
    "position_y": 0.2,
    "width": 0.03,
    "height": 0.6,
    "title_font_size": 20,
    "label_font_size": 16,
    "n_labels": 4,
    "color": "black",
    "fmt": "%.1f",
}

try:… (click to toggle)

try:
    pv.set_jupyter_backend("static")
except Exception as e:
    print("PyVista backend not set:", e)

DFN_2D = pv.read(Path(out_dir, output_mesh))
x_inlet = pv.read(Path(out_dir, "x_inlet.vtu"))
x_outlet = pv.read(Path(out_dir, "x_outlet.vtu"))


def plot_scalar_field(mesh, inlet, outlet, scalar_data, title, cmap, filename):
    scalar_bar_args = base_scalar_bar_args.copy()
    scalar_bar_args.update({"title": title, "fmt": "%.1e", "n_labels": 5})

    plotter = pv.Plotter(off_screen=True)
    plotter.add_mesh(
        mesh,
        scalars=scalar_data,
        cmap=cmap,
        show_edges=False,
        scalar_bar_args=scalar_bar_args,
    )
    plotter.add_mesh(mesh.outline(), color="black", line_width=2)
    plotter.add_mesh(inlet, color="blue", line_width=5, render_lines_as_tubes=True)
    plotter.add_mesh(outlet, color="red", line_width=5, render_lines_as_tubes=True)
    plotter.show_axes()
    plotter.enable_parallel_projection()
    plotter.view_isometric()
    plotter.screenshot(Path(out_dir, filename))
    plotter.show()
    plotter.close()


def plot_isotropic_permeability(
    mesh,
    inlet,
    outlet,
    out_dir,
    cmap="viridis",
    image_name="permeability_isotropic.png",
):
    k_iso = mesh.cell_data["permeability_ic"]
    k_log = np.log10(np.clip(k_iso, 1e-20, None))
    mesh.cell_data["log_k_iso"] = k_log

    scalar_bar_args = base_scalar_bar_args.copy()
    scalar_bar_args["title"] = "log₁₀(Permeability) [log₁₀(m²)]"

    plotter = pv.Plotter(off_screen=True)
    plotter.add_mesh(
        mesh,
        scalars="log_k_iso",
        cmap=cmap,
        show_edges=False,
        scalar_bar_args=scalar_bar_args,
    )
    plotter.add_mesh(mesh.outline(), color="black", line_width=2)
    plotter.add_mesh(inlet, color="blue", line_width=6, render_lines_as_tubes=True)
    plotter.add_mesh(outlet, color="red", line_width=6, render_lines_as_tubes=True)
    plotter.view_isometric()
    plotter.enable_parallel_projection()
    plotter.show_axes()
    plotter.screenshot(Path(out_dir, image_name))
    plotter.show()
    plotter.close()


plot_isotropic_permeability(
    mesh=DFN_2D,
    inlet=x_inlet,
    outlet=x_outlet,
    out_dir=out_dir,
)

plot_scalar_field(
    mesh=DFN_2D,
    inlet=x_inlet,
    outlet=x_outlet,
    scalar_data="width_ic",
    title="Fracture Width (m)",
    cmap="plasma",
    filename="width_ic.png",
)

png

Project file preparation

def add_component_to_model(… (click to toggle)

def add_component_to_model(
    model, name, pore_diff="1e-9", retardation="1.0", decay="0.0"
):
    model.add_element(".//media/medium/phases/phase/components", "component")
    component_xpath = ".//media/medium/phases/phase/components/component[last()]"
    model.add_element(component_xpath, "name", name)
    model.add_element(component_xpath, "properties")

    for prop, param in {
        "pore_diffusion": f"pore_diff_{name}",
        "retardation_factor": f"retard_{name}",
        "decay_rate": f"decay_{name}",
    }.items():
        model.add_element(f"{component_xpath}/properties", "property")
        prop_xpath = f"{component_xpath}/properties/property[last()]"
        model.add_element(prop_xpath, "name", prop)
        model.add_element(prop_xpath, "type", "Parameter")
        model.add_element(prop_xpath, "parameter_name", param)

    model.add_element("./process_variables", "process_variable")
    pv_xpath = "./process_variables/process_variable[last()]"

    model.add_element(pv_xpath, "name", name)
    model.add_element(pv_xpath, "components", "1")
    model.add_element(pv_xpath, "order", "1")
    model.add_element(pv_xpath, "initial_condition", "c0")
    model.add_element(pv_xpath, "boundary_conditions")
    model.add_element(f"{pv_xpath}/boundary_conditions", "boundary_condition")
    bc_xpath = f"{pv_xpath}/boundary_conditions/boundary_condition[last()]"
    model.add_element(bc_xpath, "type", "Dirichlet")
    mesh_name_bc = f"borehole_right_Layer{n_fractures}"
    model.add_element(bc_xpath, "mesh", mesh_name_bc)
    model.add_element(bc_xpath, "parameter", "c_bottom")

    model.add_element(".//processes/process/process_variables", "concentration", name)

    for pname, val in {
        f"pore_diff_{name}": pore_diff,
        f"retard_{name}": retardation,
        f"decay_{name}": decay,
    }.items():
        model.add_element(".//parameters", "parameter")
        p_xpath = ".//parameters/parameter[last()]"
        model.add_element(p_xpath, "name", pname)
        model.add_element(p_xpath, "type", "Constant")
        model.add_element(p_xpath, "value", val)


def update_reltols_for_components(model, total_components, tol="1e-12"):
    reltols_xpath = ".//time_loop/processes/process/convergence_criterion/reltols"
    reltol_values = " ".join([tol] * total_components)
    model.replace_text(reltol_values, xpath=reltols_xpath)


def update_project_parameters(project: ot.Project, params: dict):
    update_map = {
        "prefix": ".//time_loop/output/prefix",
        "initial_pressure_expression": ".//parameters/parameter[name='p0']/expression",
        "inlet_pressure_expression": ".//parameters/parameter[name='pinlet']/expression",
        "outlet_pressure_expression": ".//parameters/parameter[name='poutlet']/expression",
        "inlet_concentration_value": ".//parameters/parameter[name='c_bottom']/value",
        "porosity_value": ".//parameters/parameter[name='porosity_parameter']/value",
        "fracture_thickness_value": ".//parameters/parameter[name='fracture_thickness_const']/value",
        "decay_Si_value": ".//parameters/parameter[name='decay']/value",
        "t_end": ".//time_loop/processes/process/time_stepping/t_end",
        "maximum_dt": ".//time_loop/processes/process/time_stepping/maximum_dt",
    }

    for key, xpath in update_map.items():
        if key in params:
            try:
                project.replace_text(params[key], xpath=xpath)
                display(
                    Markdown(
                        f"**Success:** Parameter `{key}` updated to **{params[key]}**"
                    )
                )
            except Exception as e:
                msg = f"Failed to update parameter `{key}`: {e}"
                display(Markdown(f"**Error:** {msg}"))
                raise RuntimeError(msg) from e

    if "fracture_permeability_values" in params:
        try:
            project.replace_text(
                params["fracture_permeability_values"].strip(),
                xpath=".//parameters/parameter[name='kappa1_frac']/values",
            )
            display(
                Markdown(
                    "**Success:** Fracture permeability values updated successfully."
                )
            )
        except Exception as e:
            msg = f"Failed to update fracture permeability: {e}"
            display(Markdown(f"**Error:** {msg}"))
            raise RuntimeError(msg) from e

    main_proc_xpath = "processes/process[1]"
    project.remove_element(f"{main_proc_xpath}/numerical_stabilization")

    if "numerical_stabilization" in params:
        stab = params["numerical_stabilization"]
        stype = stab.get("type")
        if not stype:
            error_msg = "'type' is required in numerical_stabilization"
            raise ValueError(error_msg)

        project.add_element(main_proc_xpath, "numerical_stabilization")
        ns_xpath = f"{main_proc_xpath}/numerical_stabilization[last()]"
        project.add_element(ns_xpath, "type", stype)

        if stype == "FullUpwind":
            cutoff = stab.get("cutoff_velocity", "0.0")
            project.add_element(ns_xpath, "cutoff_velocity", cutoff)

        elif stype == "IsotropicDiffusion":
            tuning = stab.get("tuning_parameter")
            cutoff = stab.get("cutoff_velocity")
            if tuning is None or cutoff is None:
                error_msg = "'tuning_parameter' and 'cutoff_velocity' required for IsotropicDiffusion"
                raise ValueError(error_msg)
            project.add_element(ns_xpath, "tuning_parameter", tuning)
            project.add_element(ns_xpath, "cutoff_velocity", cutoff)

        elif stype == "FluxCorrectedTransport":
            pass

        else:
            error_msg = f"Unsupported stabilization type: {stype}"
            raise ValueError(error_msg)

        display(
            Markdown(
                f"**Success:** Configured numerical stabilization in main processes with type `{stype}`"
            )
        )

    n = params["n_fractures"]
    bc_mesh_xpath = ".//process_variables/process_variable[name='Si']/boundary_conditions/boundary_condition/mesh"
    project.replace_text(f"borehole_right_Layer{n}", xpath=bc_mesh_xpath)

    display(
        Markdown(f"**Success:** Updated BC mesh for `Si` to `borehole_right_Layer{n}`")
    )

    left_mesh = f"borehole_left_Layer{n + 1}.vtu"
    right_mesh = f"borehole_right_Layer{n}.vtu"

    mesh_base_xpath = ".//meshes/mesh"
    project.replace_text(left_mesh, xpath=f"{mesh_base_xpath}[last()-1]")
    project.replace_text(right_mesh, xpath=f"{mesh_base_xpath}[last()]")
    display(
        Markdown(f"**Success:** Updated borehole meshes: `{left_mesh}`, `{right_mesh}`")
    )

user_parameters = {… (click to toggle)

user_parameters = {
    "prefix": "DFN_HC",
    "initial_pressure_expression": "1000*9.81*(500-z)",
    "inlet_pressure_expression": "1000*9.81*(500-z)+2.943e5",
    "outlet_pressure_expression": "1000*9.81*(500-z)",
    "inlet_concentration_value": "1",
    "porosity_value": "0.05",
    "decay_Si_value": "1e-12",
    "t_end": "1e11",
    "maximum_dt": "5e9",
    "numerical_stabilization": {
        "type": "FluxCorrectedTransport",
    },
    "n_fractures": n_fractures,
}

project_file = Path("DFN_HC.prj")… (click to toggle)

project_file = Path("DFN_HC.prj")

project = ot.Project(
    input_file=project_file, output_file=Path(f"{out_dir}/DFN_HC_final.prj")
)
add_component_to_model(
    project, "Cs", pore_diff="1e-9", retardation="1.", decay="1.e-10"
)
update_reltols_for_components(
    project, total_components=3
)  # total = 1 pressure + 2 components
update_project_parameters(project, user_parameters)
project.write_input()

Success: Parameter prefix updated to DFN_HC

Success: Parameter initial_pressure_expression updated to 10009.81(500-z)

Success: Parameter inlet_pressure_expression updated to 10009.81(500-z)+2.943e5

Success: Parameter outlet_pressure_expression updated to 10009.81(500-z)

Success: Parameter inlet_concentration_value updated to 1

Success: Parameter porosity_value updated to 0.05

Success: Parameter decay_Si_value updated to 1e-12

Success: Parameter t_end updated to 1e11

Success: Parameter maximum_dt updated to 5e9

Success: Configured numerical stabilization in main processes with type FluxCorrectedTransport

Success: Updated BC mesh for Si to borehole_right_Layer10

Success: Updated borehole meshes: borehole_left_Layer11.vtu, borehole_right_Layer10.vtu

Run simulations

project.run_model(args=f"-o {out_dir} -m {out_dir}", logfile=Path(out_dir, "run.log"))

 ogs -o /var/lib/gitlab-runner/builds/t1_R7hwLX/1/ogs/build/release-all/Tests/Data/Parabolic/ComponentTransport/DFN_PorePy/DFNbyPorePy_to_OGS -m /var/lib/gitlab-runner/builds/t1_R7hwLX/1/ogs/build/release-all/Tests/Data/Parabolic/ComponentTransport/DFN_PorePy/DFNbyPorePy_to_OGS /var/lib/gitlab-runner/builds/t1_R7hwLX/1/ogs/build/release-all/Tests/Data/Parabolic/ComponentTransport/DFN_PorePy/DFNbyPorePy_to_OGS/DFN_HC_final.prj
['ogs', '-o', '/var/lib/gitlab-runner/builds/t1_R7hwLX/1/ogs/build/release-all/Tests/Data/Parabolic/ComponentTransport/DFN_PorePy/DFNbyPorePy_to_OGS', '-m', '/var/lib/gitlab-runner/builds/t1_R7hwLX/1/ogs/build/release-all/Tests/Data/Parabolic/ComponentTransport/DFN_PorePy/DFNbyPorePy_to_OGS', '/var/lib/gitlab-runner/builds/t1_R7hwLX/1/ogs/build/release-all/Tests/Data/Parabolic/ComponentTransport/DFN_PorePy/DFNbyPorePy_to_OGS/DFN_HC_final.prj']

<Popen: returncode: 0 args: ' ogs -o /var/lib/gitlab-runner/builds/t1_R7hwLX...>

Post-processing

ms = ot.MeshSeries(f'{out_dir}/{user_parameters["prefix"]}.pvd')
mesh = ms[-1]

plotter = pv.Plotter(off_screen=True)… (click to toggle)

plotter = pv.Plotter(off_screen=True)
scalar_bar_args = base_scalar_bar_args.copy()
scalar_bar_args["title"] = "Pressure [Pa]"
plotter.add_mesh(
    mesh,
    scalars="pressure",
    cmap="Blues",
    show_edges=True,
    opacity=1.0,
    scalar_bar_args=scalar_bar_args,
)


plotter.add_mesh(mesh.outline(), color="black", line_width=2)
plotter.show_axes()
plotter.enable_parallel_projection()
plotter.view_isometric()

output_path = out_dir.joinpath("pressure.png")
plotter.screenshot(str(output_path))
plotter.show()
plotter.close()

png

magnitude = np.linalg.norm(mesh["darcy_velocity"], axis=1)… (click to toggle)

magnitude = np.linalg.norm(mesh["darcy_velocity"], axis=1)
magnitude[magnitude <= 0] = 1e-12
log_magnitude = np.log10(magnitude)
mesh["darcy_velocity_magnitude"] = magnitude
mesh["log_darcy_velocity_magnitude"] = log_magnitude
mesh.set_active_vectors("darcy_velocity")

vmax = magnitude.max()
domain_size = max(mesh.bounds[1::2]) - min(mesh.bounds[::2])
scale_factor = domain_size * 0.1 / vmax if vmax > 0 else 1.0

arrows = mesh.glyph(
    scale="darcy_velocity_magnitude", orient="darcy_velocity", factor=scale_factor
)

scalar_bar_args = base_scalar_bar_args.copy()
scalar_bar_args["title"] = "log₁₀(Darcy Velocity) [m/s]"

plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(
    mesh,
    scalars="log_darcy_velocity_magnitude",
    cmap="viridis",
    opacity=1.0,
    show_edges=False,
    scalar_bar_args=scalar_bar_args,
)
plotter.add_mesh(arrows, color="black", label="Velocity Vectors")
plotter.add_mesh(mesh.outline(), color="black", line_width=2)
plotter.show_axes()
plotter.enable_parallel_projection()
plotter.view_isometric()
output_path = out_dir.joinpath("full_mesh_velocity_log_arrows.png")
plotter.screenshot(str(output_path))
plotter.show()
plotter.close()

png

plotter = pv.Plotter(off_screen=True)… (click to toggle)

plotter = pv.Plotter(off_screen=True)
scalar_bar_args = base_scalar_bar_args.copy()
scalar_bar_args["title"] = "Si Concentration [mol/m³]"

plotter.add_mesh(
    mesh,
    scalars="Si",
    cmap="plasma",
    clim=(0, 1),
    show_edges=False,
    opacity=1.0,
    scalar_bar_args=scalar_bar_args,
)

plotter.add_mesh(mesh.outline(), color="black", line_width=2)
plotter.show_axes()
plotter.enable_parallel_projection()
plotter.view_isometric()

output_path = out_dir.joinpath("concentration.png")
plotter.screenshot(str(output_path))
plotter.show()
plotter.close()

png

plotter = pv.Plotter(off_screen=True)… (click to toggle)

plotter = pv.Plotter(off_screen=True)
scalar_bar_args = base_scalar_bar_args.copy()
scalar_bar_args["title"] = "Cs Concentration [mol/m³]"

plotter.add_mesh(
    mesh,
    scalars="Cs",
    clim=(0, 1),
    cmap="cividis",
    show_edges=False,
    opacity=1.0,
    scalar_bar_args=scalar_bar_args,
)

plotter.add_mesh(mesh.outline(), color="black", line_width=2)
plotter.show_axes()
plotter.enable_parallel_projection()
plotter.view_isometric()

output_path = out_dir.joinpath("concentration_Cs.png")
plotter.screenshot(str(output_path))
plotter.show()
plotter.close()

png

Breakthrough curve (BTC)

the temporal evolution of tracer concentration at a downstream monitoring point in a flow-through porous medium,

mesh_series = ot.MeshSeries(f'{out_dir}/{user_parameters["prefix"]}.pvd')… (click to toggle)

mesh_series = ot.MeshSeries(f'{out_dir}/{user_parameters["prefix"]}.pvd')
print("Scalar fields available in mesh:", mesh.array_names)

observation_points = np.array(
    [
        [80, 50, 90],
        [80, 50, 80],
    ]
)

labels = [f"Point {i}: {pt}" for i, pt in enumerate(observation_points)]


ms_days = mesh_series.scale(time="d")
var = ot.variables.Scalar("Si")

ms_probe = ms_days.probe(points=observation_points)
fig_si = ms_probe.plot_line("time", var, labels=labels)
fig_si.axes[0].set_xscale("log")

fig_cs = ms_probe.plot_line("time", "Cs", labels=labels)
fig_cs.axes[0].set_xscale("log")

Scalar fields available in mesh: ['Cs', 'OGS_VERSION', 'CsFlowRate', 'LiquidMassFlowRate', 'Si', 'SiFlowRate', 'darcy_velocity', 'pressure', 'darcy_velocity_magnitude', 'log_darcy_velocity_magnitude', 'MaterialIDs', 'permeability_ic', 'porosity_avg', 'velocity', 'width_ic']

png

def _ensure_point_field(m, field_name):… (click to toggle)

def _ensure_point_field(m, field_name):
    if field_name in m.point_data:
        return m
    if field_name in m.cell_data:
        return m.cell_data_to_point_data()
    msg = f"{field_name} not found in point_data or cell_data"
    raise KeyError(msg)


def _nn_fill(masked_vals, valid_mask, tgt_pts, src_pts, src_vals):
    nan_idx = ~valid_mask
    if not np.any(nan_idx):
        return masked_vals
    tree = cKDTree(src_pts)
    _, j = tree.query(tgt_pts[nan_idx], k=1)
    masked_vals[nan_idx] = src_vals[j]
    return masked_vals


last_mesh = mesh_series[-1]
pressure_last = np.asarray(last_mesh.point_data["pressure"])

ref_mesh_path = Path("DFN_HC_ts_51_t_100000000000.000000.vtu")
if not ref_mesh_path.exists():
    available = list(Path().glob("DFN_HC_ts_*.vtu"))
    msg = (
        f"Reference mesh not found: {ref_mesh_path.resolve()}\n"
        f"This file is generated by DFNbyPorePy.py during the test setup.\n"
        f"Available files in current directory: {available if available else 'none'}"
    )
    raise FileNotFoundError(msg)

ref_mesh = pv.read(ref_mesh_path)


def _sample_ref_field_on_last_mesh(field_name):
    ref_mesh_field = _ensure_point_field(ref_mesh, field_name)
    sampled = last_mesh.sample(ref_mesh_field)
    ref_values_on_last = np.asarray(sampled.point_data[field_name])

    mask = sampled.point_data.get("vtkValidPointMask", None)
    if mask is not None:
        mask = np.asarray(mask).astype(bool)
        if not mask.all():
            ref_values_on_last = _nn_fill(
                ref_values_on_last,
                mask,
                last_mesh.points,
                ref_mesh_field.points,
                np.asarray(ref_mesh_field.point_data[field_name]),
            )
    return ref_values_on_last


pressure_ref_on_last = _sample_ref_field_on_last_mesh("pressure")

np.testing.assert_allclose(
    actual=pressure_last,
    desired=pressure_ref_on_last,
    rtol=5e-6,
    atol=1e-8,
    equal_nan=False,
)
print("Pressure fields match")

for field_name in ["Si", "Cs"]:
    if field_name in last_mesh.point_data or field_name in last_mesh.cell_data:
        last_values = np.asarray(
            _ensure_point_field(last_mesh, field_name).point_data[field_name]
        )
        ref_values = _sample_ref_field_on_last_mesh(field_name)
        np.testing.assert_allclose(
            actual=last_values,
            desired=ref_values,
            rtol=1e-3,
            atol=1e-8,
            equal_nan=False,
        )
        print(f"{field_name} fields match")

Pressure fields match
Si fields match
Cs fields match

References

PorePy – An open-source simulation tool for fractured and porous media. Available at: https://github.com/pmgbergen/porepy
Keilegavlen, E., Berre, I., Fumagalli, A., Stefansson, I., and Edwards, M.G. (2021). PorePy: An open-source simulation tool for flow and transport in deformable fractured rocks. Computational Geosciences, 25, 165–187.
Buchwald, J., Kolditz, O., & Nagel, T. (2021). ogs6py and VTUinterface: streamlining OpenGeoSys workflows in Python. Journal of Open Source Software, 6(67), 3673.
Cvetkovic, V., & Frampton, A. (2010). Transport and retention from single to multiple fractures in crystalline rock at Äspö (Sweden): Fracture network simulations and generic retention model. Water Resources Research, 46, W05506.
Watanabe, N., & Kolditz, O. (2015). Numerical stability analysis of two-dimensional solute transport along a discrete fracture in a porous rock matrix. Water Resources Research, 51(7), 5855-5868.
Chen, C., Binder, M., Oppelt, L., Hu, Y., Engelmann, C., Arab, A., Xu, W., Scheytt, T., & Nagel, T. (2024). Modeling of heat and solute transport in a fracture-matrix mine thermal energy storage system and energy storage performance evaluation. Journal of Hydrology, 636(May), 131335.
Chen, R., Xu, W., Chen, Y., Nagel, T., Chen, C., Hu, Y., Li, J., & Zhuang, D. (2024). Influence of heterogeneity on dissolved CO2 migration in a fractured reservoir. Environmental Earth Sciences, 83(16), 463.
Chaudhry, A. A., Zhang, C., Ernst, O. G., & Nagel, T. (2025). Effects of inhomogeneity and statistical and material anisotropy on THM simulations. Reliability Engineering & System Safety, 110921.
Zheng, Z., Wang, X., Flügge, J., & Nagel, T. (2025). A stochastic modeling framework for radionuclide migration from deep geological repositories considering spatial variability. Advances in Water Resources, 203, 105003.

Reactive transport and retention in stochastically generated DFNs: A PorePy–OpenGeoSys Workflow

Reactive transport and retention in stochastically generated DFNs: A PorePy–OpenGeoSys Workflow

Why is solute retention critical?

How to study the solute retention?

How to model the transport in fractured rock?

How to numerically model it?

Problem description

OpenGeosys

Preparing the mesh (.vtu) for OpenGeosys simulation

Boundary and materials extraction using OpenGeoSys tools

Plot DFN and boundaries

Generate and assign random fracture properties with Gaussian distributions

Project file preparation

Run simulations

Post-processing

Breakthrough curve (BTC)

References

Preparing the mesh (`.vtu`) for OpenGeosys simulation