check_profugus.py 3.59 KB
Newer Older
Hamilton, Steven P's avatar
Hamilton, Steven P committed
1
2
3
4
5
6

import h5py
import numpy as np
import os.path

def check_solution(path_to_results):
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
    # Load reference flux
    ref_flux_filename = os.path.join(path_to_results,'c5g7_3d_flux_ref.h5')
    assert os.path.isfile(ref_flux_filename)
    ref_flux_file = h5py.File(ref_flux_filename,'r')
    ref_flux     = np.array(ref_flux_file['flux_mean'])
    ref_flux_std = np.array(ref_flux_file['flux_std_dev'])

    # Load reference k
    ref_k_filename = os.path.join(path_to_results,'c5g7_3d_output_ref.h5')
    assert os.path.isfile(ref_k_filename)
    ref_k_file = h5py.File(ref_k_filename,'r')
    ref_k_dset     = ref_k_file['/keff/mean']
    ref_k_var_dset = ref_k_file['/keff/variance']
    ref_k = ref_k_dset[0]
    ref_k_std = np.sqrt(ref_k_var_dset[0])

    print('\nReference k = {0:.6f} +/- {1:.3e}'.format(ref_k,ref_k_std))

    # Load flux file for comparison
    flux_filename = os.path.join(path_to_results,'c5g7_3d_flux.h5')
    assert os.path.isfile(flux_filename)
    flux_file = h5py.File(flux_filename,'r')
    flux     = np.array(flux_file['flux_mean'])
    flux_std = np.array(flux_file['flux_std_dev'])

    # Load k file for comparison
    k_filename = os.path.join(path_to_results,'c5g7_3d_output.h5')
    assert os.path.isfile(k_filename)
    k_file = h5py.File(k_filename,'r')
    k_dset     = k_file['/keff/mean']
    k_var_dset = k_file['/keff/variance']
    k = k_dset[0]
    k_std = np.sqrt(k_var_dset[0])

    print('New k = {0:.6f} +/- {1:.3e}'.format(k,k_std))

    passes = True

    # Test to 4 standard deviations, should *almost* never fail
    safety_factor = 4

    # Check k values
    kdiff = np.abs(k - ref_k)
    tol = safety_factor*np.sqrt(k_std*k_std + ref_k_std*ref_k_std)
    print('Difference in k = {0:.2e} with tolerance of {1:.2e}'.format(kdiff,tol))
    passes = kdiff < tol

    if passes:
        'Eigenvalue difference is acceptable'
    else:
        'Eigenvalue difference is outside expected tolerance'

    # Cycle-to-cycle correlations cause reported standard deviation to
    # be under-reported for tally values.  Increase safety factor on flux
    # checks to eliminate false positives.
    safety_factor = 10

    # Check flux
    N = len(ref_flux)
    assert N == len(ref_flux_std)
    assert N == len(flux)
    assert N == len(flux_std)

    failed_indices = []
    max_ratio = 0.0
    for i in xrange(N):
        ref_val = ref_flux[i]
        ref_std = ref_flux_std[i]
        val     = flux[i]
        std     = flux_std[i]
        tol = safety_factor*np.sqrt(ref_std*ref_std + std*std)
        flux_diff = np.abs(val - ref_val)
        ratio = flux_diff / tol
        max_ratio = max(ratio,max_ratio)
        if flux_diff > tol:
            failed_indices.append(i)

    num_failed = len(failed_indices)

    print('\nMaximum relative error is {0:.3f} times tolerance'.format(max_ratio))

    # Allow one outlier
    if num_failed > 1:
        passes = False

    if num_failed > 0:
        print('\nThe following flux values were outside of expected tolerance:')
        print(' Index   Reference   Value       Difference Tolerance')
        for ind in failed_indices:
            ref_val = ref_flux[ind]
            ref_std = ref_flux_std[ind]
            val     = flux[ind]
            std     = flux_std[ind]
            tol = safety_factor*np.sqrt(ref_std*ref_std + std*std)
            flux_diff = np.abs(val - ref_val)
            print('{0:5}    {1:.4e}  {2:.4e}  {3:.3e}  {4:.3e}'.format(
                ind,ref_val,val,flux_diff,tol))

    # Print pass/fail status
    print('\n'),
    if passes:
        print('Solution PASSED')
        return 0
    else:
        print('Solution FAILED')
        return 1
Hamilton, Steven P's avatar
Hamilton, Steven P committed
113