-
Notifications
You must be signed in to change notification settings - Fork 21
Expand file tree
/
Copy pathplotmodel.py
More file actions
1782 lines (1499 loc) · 65.1 KB
/
plotmodel.py
File metadata and controls
1782 lines (1499 loc) · 65.1 KB
1
2
3
4
5
6
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
from __future__ import annotations
from ast import literal_eval
from collections import defaultdict
import copy
from ctypes import c_int32, c_char_p
from dataclasses import dataclass
import hashlib
import itertools
import pickle
from typing import Literal, Tuple, Optional, Dict
from PySide6.QtWidgets import QItemDelegate, QColorDialog, QLineEdit, QMessageBox
from PySide6.QtCore import (QAbstractTableModel, QModelIndex, Qt, QSize, QEvent,
QObject, Signal, Slot, QThread, QEventLoop, QTimer)
from PySide6.QtGui import QColor
import openmc
import openmc.lib
import numpy as np
from . import __version__
from .statepointmodel import StatePointModel
from .plot_colors import random_rgb
ID, NAME, COLOR, COLORLABEL, MASK, HIGHLIGHT = range(6)
_VOID_REGION = -1
_NOT_FOUND = -2
_OVERLAP = -3
_MODEL_PROPERTIES = ('temperature', 'density')
_PROPERTY_INDICES = {'temperature': 0, 'density': 1}
_REACTION_UNITS = 'reactions/source'
_PRODUCTION_UNITS = 'particles/source'
_ENERGY_UNITS = 'eV/source'
_REACTION_UNITS_VOL = 'reactions/cm³/source'
_PRODUCTION_UNITS_VOL = 'particles/cm³/source'
_ENERGY_UNITS_VOL = 'eV/cm³/source'
_SPATIAL_FILTERS = (openmc.UniverseFilter,
openmc.MaterialFilter,
openmc.CellFilter,
openmc.DistribcellFilter,
openmc.CellInstanceFilter,
openmc.MeshFilter,
openmc.MeshMaterialFilter)
_PRODUCTIONS = ('delayed-nu-fission', 'prompt-nu-fission', 'nu-fission',
'nu-scatter', 'H1-production', 'H2-production',
'H3-production', 'He3-production', 'He4-production')
_ENERGY_SCORES = {'heating', 'heating-local', 'kappa-fission',
'fission-q-prompt', 'fission-q-recoverable',
'damage-energy'}
_SCORE_UNITS = {p: _PRODUCTION_UNITS for p in _PRODUCTIONS}
_SCORE_UNITS['flux'] = 'particle-cm/source'
_SCORE_UNITS['current'] = 'particle/source'
_SCORE_UNITS['events'] = 'events/source'
_SCORE_UNITS['inverse-velocity'] = 'particle-s/source'
_SCORE_UNITS['decay-rate'] = 'particle/s/source'
_SCORE_UNITS.update({s: _ENERGY_UNITS for s in _ENERGY_SCORES})
_SCORE_UNITS_VOL = {p: _PRODUCTION_UNITS_VOL for p in _PRODUCTIONS}
_SCORE_UNITS_VOL['flux'] = 'particle/cm²/source'
_SCORE_UNITS_VOL['current'] = 'particle/cm³/source'
_SCORE_UNITS_VOL['events'] = 'events/cm³/source'
_SCORE_UNITS_VOL['inverse-velocity'] = 'particle-s/cm³/source'
_SCORE_UNITS_VOL['decay-rate'] = 'particle/s/cm³/source'
_SCORE_UNITS.update({s: _ENERGY_UNITS_VOL for s in _ENERGY_SCORES})
_TALLY_VALUES = {'Mean': 'mean',
'Std. Dev.': 'std_dev',
'Rel. Error': 'rel_err'}
TallyValueType = Literal['mean', 'std_dev', 'rel_err']
def hash_file(path):
# return the md5 hash of a file
h = hashlib.md5()
with path.open('rb') as file:
chunk = 0
while chunk != b'':
# read 32768 bytes at a time
chunk = file.read(32768)
h.update(chunk)
return h.hexdigest()
def hash_model(model_path):
"""Get hash values for materials.xml and geometry.xml (or model.xml)"""
if model_path.is_file():
mat_xml_hash = hash_file(model_path)
geom_xml_hash = ""
elif (model_path / 'model.xml').exists():
mat_xml_hash = hash_file(model_path / 'model.xml')
geom_xml_hash = ""
else:
mat_xml_hash = hash_file(model_path / 'materials.xml')
geom_xml_hash = hash_file(model_path / 'geometry.xml')
return mat_xml_hash, geom_xml_hash
@dataclass(frozen=True)
class PlotRequest:
view_snapshot: "PlotView"
view_params: Dict[str, object]
@dataclass(frozen=True)
class PlotWorkItem:
view_params: Dict[str, object]
class PlotWorker(QObject):
finished = Signal(object, object, object)
error = Signal(str)
@Slot(object)
def generate_maps(self, work_item: PlotWorkItem):
try:
params = work_item.view_params
# Determine if we need filter bins for MeshMaterialFilter tally
filter_cpp = None
if params.get("filter_id") is not None:
filter_cpp = openmc.lib.filters[params["filter_id"]]
# Single call replaces id_map + property_map + get_plot_bins
geom_data, property_data = openmc.lib.slice_plot(
origin=params["origin"],
width=(params["width"], params["height"]),
basis=params["basis"],
pixels=(params["h_res"], params["v_res"]),
color_overlaps=params["color_overlaps"],
level=params["level"],
filter=filter_cpp,
)
self.finished.emit(work_item.view_params, geom_data, property_data)
except Exception as exc:
self.error.emit(str(exc))
class PlotManager(QObject):
plot_started = Signal()
plot_queued = Signal()
plot_finished = Signal(object, object, object, object)
plot_error = Signal(str)
plot_idle = Signal()
work_requested = Signal(object)
def __init__(self):
super().__init__()
self._thread = QThread()
self._worker = PlotWorker()
self._worker.moveToThread(self._thread)
self._thread.finished.connect(self._worker.deleteLater)
self.work_requested.connect(self._worker.generate_maps)
self._worker.finished.connect(self._on_worker_finished)
self._worker.error.connect(self._on_worker_error)
self._thread.start()
self._pending_request = None
self._in_flight_request = None
self._latest_view_params = None
@property
def latest_view_params(self):
return self._latest_view_params
@property
def is_busy(self):
return self._in_flight_request is not None
@property
def has_pending(self):
return self._pending_request is not None
def enqueue(self, view_snapshot, view_params):
request = PlotRequest(view_snapshot, view_params)
self._latest_view_params = view_params
if self._in_flight_request is None:
self._pending_request = request
self._start_next()
return
self._pending_request = request
self.plot_queued.emit()
def set_latest_view_params(self, view_params):
self._latest_view_params = view_params
def clear_pending(self):
self._pending_request = None
def wait_for_idle(self, timeout_ms: Optional[int] = None):
if not self.is_busy and not self.has_pending:
return True
loop = QEventLoop()
timer = None
if timeout_ms is not None:
timer = QTimer()
timer.setSingleShot(True)
timer.timeout.connect(loop.quit)
def on_idle():
loop.quit()
self.plot_idle.connect(on_idle)
if timer is not None:
timer.start(timeout_ms)
loop.exec()
self.plot_idle.disconnect(on_idle)
if timer is not None:
timer.stop()
return not self.is_busy and not self.has_pending
def shutdown(self):
self._pending_request = None
self._in_flight_request = None
if self._thread.isRunning():
self._thread.quit()
self._thread.wait()
def _start_next(self):
if self._pending_request is None:
return
self._in_flight_request = self._pending_request
self._pending_request = None
self.plot_started.emit()
work_item = PlotWorkItem(self._in_flight_request.view_params)
self.work_requested.emit(work_item)
@Slot(object, object, object)
def _on_worker_finished(self, view_params, geom_data, property_data):
request = self._in_flight_request
self._in_flight_request = None
if request is not None:
self.plot_finished.emit(request.view_snapshot,
view_params, geom_data, property_data)
if self._pending_request is not None:
self._start_next()
else:
self.plot_idle.emit()
@Slot(str)
def _on_worker_error(self, message):
self._in_flight_request = None
self.plot_error.emit(message)
if self._pending_request is not None:
self._start_next()
else:
self.plot_idle.emit()
class PlotModel:
"""Geometry and plot settings for OpenMC Plot Explorer model
Parameters
----------
use_settings_pkl : bool
If True, use plot_settings.pkl file to reload settings
model_path : pathlib.Path
Path to model XML file or directory
default_res : int
Default resolution as the number of pixels in each direction
Attributes
----------
geom : openmc.Geometry
OpenMC Geometry of the model
modelCells : collections.OrderedDict
Dictionary mapping cell IDs to openmc.Cell instances
modelMaterials : collections.OrderedDict
Dictionary mapping material IDs to openmc.Material instances
ids : NumPy int array (v_res, h_res, 1)
Mapping of plot coordinates to cell/material ID by pixel
geom_data : NumPy int32 array (v_res, h_res, 3) or (v_res, h_res, 4)
Geometry data with cell IDs, instances, material IDs, and optionally filter bins
property_data : Numpy float array (v_res, h_res, 2)
Property data with cell temperatures and material densities
image : NumPy int array (v_res, h_res, 3)
The current RGB image data
statepoint : StatePointModel
Simulation data model used to display tally results
applied_filters : tuple of ints
IDs of the applied filters for the displayed tally
previousViews : list of PlotView instances
List of previously created plot view settings used to undo
changes made in plot explorer
subsequentViews : list of PlotView instances
List of undone plot view settings used to redo changes made
in plot explorer
sourceSitesTolerance : float
Tolerance for source site plotting (default 0.1 cm)
sourceSitesColor : tuple of 3 int
RGB color for source site plotting (default blue)
sourceSitesVisible : bool
Whether to plot source sites (default True)
sourceSites : Source sites to plot
Set of source locations to plot
defaultView : PlotView
Default settings for given geometry
currentView : PlotView
Currently displayed plot settings in plot explorer
activeView : PlotView
Active state of settings in plot explorer, which may or may not
have unapplied changes
"""
def __init__(self, use_settings_pkl, model_path, default_res):
""" Initialize PlotModel class attributes """
# Retrieve OpenMC Cells/Materials
self.modelCells = openmc.lib.cells
self.modelMaterials = openmc.lib.materials
self.max_universe_levels = openmc.lib._coord_levels()
# Cell/Material ID by coordinates
self.ids = None
# Return values from slice_plot
self.geom_data = None
self.property_data = None
self.map_view_params = None
self.version = __version__
# default statepoint value
self._statepoint = None
# default tally/filter info
self.appliedFilters = ()
self.appliedScores = ()
self.appliedNuclides = ()
self.previousViews = []
self.subsequentViews = []
self.defaultView = self.getDefaultView(default_res)
# Source site defaults
self.sourceSitesApplyTolerance = False
self.sourceSitesTolerance = 0.1 # cm
self.sourceSitesColor = (0, 0, 255)
self.sourceSitesVisible = True
self.sourceSites = None
if model_path.is_file():
settings_pkl = model_path.with_name('plot_settings.pkl')
else:
settings_pkl = model_path / 'plot_settings.pkl'
if use_settings_pkl and settings_pkl.is_file():
with settings_pkl.open('rb') as file:
try:
data = pickle.load(file)
except AttributeError:
msg_box = QMessageBox()
msg = "WARNING: previous plot settings are in an incompatible format. " +\
"They will be ignored."
msg_box.setText(msg)
msg_box.setIcon(QMessageBox.Warning)
msg_box.setStandardButtons(QMessageBox.Ok)
msg_box.exec()
self.currentView = copy.deepcopy(self.defaultView)
else:
restore_domains = False
# check GUI version
if data['version'] != self.version:
print("WARNING: previous plot settings are for a different "
"version of the GUI. They will be ignored.")
wrn_msg = "Existing version: {}, Current GUI version: {}"
print(wrn_msg.format(data['version'], self.version))
view = None
else:
view = data['currentView']
# get materials.xml and geometry.xml hashes to
# restore additional settings if possible
mat_xml_hash, geom_xml_hash = hash_model(model_path)
if mat_xml_hash == data['mat_xml_hash'] and \
geom_xml_hash == data['geom_xml_hash']:
restore_domains = True
# restore statepoint file
try:
self.statepoint = data['statepoint']
except OSError:
msg_box = QMessageBox()
msg = "Could not open statepoint file: \n\n {} \n"
msg_box.setText(msg.format(self.model.statepoint.filename))
msg_box.setIcon(QMessageBox.Warning)
msg_box.setStandardButtons(QMessageBox.Ok)
msg_box.exec()
self.statepoint = None
self.currentView = PlotView(restore_view=view,
restore_domains=restore_domains,
default_res=default_res)
else:
self.currentView = copy.deepcopy(self.defaultView)
self.activeView = copy.deepcopy(self.currentView)
self.plot_manager = PlotManager()
def openStatePoint(self, filename):
self.statepoint = StatePointModel(filename, open_file=True)
@property
def statepoint(self):
return self._statepoint
@statepoint.setter
def statepoint(self, statepoint):
if statepoint is None:
self._statepoint = None
elif isinstance(statepoint, StatePointModel):
self._statepoint = statepoint
elif isinstance(statepoint, str):
self._statepoint = StatePointModel(statepoint, open_file=True)
else:
raise TypeError("Invalid statepoint object")
if self._statepoint and not self._statepoint.is_open:
self._statepoint.open()
def getDefaultView(self, default_res):
""" Generates default PlotView instance for OpenMC geometry
Centers plot view origin in every dimension if possible. Defaults
to xy basis, with height and width to accomodate full size of
geometry. Defaults to (0, 0, 0) origin with width and heigth of
25 if geometry bounding box cannot be generated.
Parameters
----------
default_res : int
Default resolution as the number of pixels in each direction
Returns
-------
default : PlotView instance
PlotView instance with default view settings
"""
lower_left, upper_right = openmc.lib.global_bounding_box()
# Check for valid bounding_box dimensions
if -np.inf not in lower_left[:2] and np.inf not in upper_right[:2]:
xcenter = (upper_right[0] + lower_left[0])/2
width = abs(upper_right[0] - lower_left[0]) * 1.005
ycenter = (upper_right[1] + lower_left[1])/2
height = abs(upper_right[1] - lower_left[1]) * 1.005
else:
xcenter, ycenter, width, height = (0.00, 0.00, 25, 25)
if lower_left[2] != -np.inf and upper_right[2] != np.inf:
zcenter = (upper_right[2] + lower_left[2])/2
else:
zcenter = 0.00
default = PlotView([xcenter, ycenter, zcenter], width, height,
default_res=default_res)
return default
def resetColors(self):
""" Reset colors to those generated in the default view """
self.activeView.cells = self.defaultView.cells
self.activeView.materials = self.defaultView.materials
def view_params_payload(self, view: "PlotView"):
vp = view.view_params
return {
"origin": tuple(vp.origin),
"width": float(vp.width),
"height": float(vp.height),
"h_res": int(vp.h_res),
"v_res": int(vp.v_res),
"basis": str(vp.basis),
"level": int(vp.level),
"color_overlaps": bool(vp.color_overlaps),
"filter_id": self.get_active_mesh_material_filter_id(view),
}
def get_active_mesh_material_filter_id(self, view: "PlotView") -> Optional[int]:
"""Return the filter ID if displaying a MeshMaterialFilter tally, else None."""
if self._statepoint is None:
return None
if not view.tallyDataVisible or view.selectedTally is None:
return None
tally = self._statepoint.tallies[view.selectedTally]
if tally.contains_filter(openmc.MeshMaterialFilter):
filter = tally.find_filter(openmc.MeshMaterialFilter)
return filter.id
return None
def can_reuse_maps(self, view: "PlotView"):
if self.geom_data is None or self.property_data is None:
return False
return self.map_view_params == self.view_params_payload(view)
def makePlot(self, view: Optional["PlotView"] = None,
geom_data=None, property_data=None):
"""Generate new plot image from active view settings"""
if view is None:
view = self.activeView
# update/call maps under 2 circumstances
# 1. this is the intial plot (geom_data/property_data are None)
# 2. The active (desired) view differs from the current view parameters
if geom_data is None or property_data is None:
if (self.currentView.view_params != view.view_params) or \
(self.geom_data is None) or (self.property_data is None):
# Determine if we need filter bins for MeshMaterialFilter tally
filter_cpp = None
filter_id = self.get_active_mesh_material_filter_id(view)
if filter_id is not None:
filter_cpp = openmc.lib.filters[filter_id]
self.geom_data, self.property_data = openmc.lib.slice_plot(
origin=view.origin,
width=(view.width, view.height),
basis=view.basis,
pixels=(view.h_res, view.v_res),
color_overlaps=view.color_overlaps,
level=view.level,
filter=filter_cpp,
)
self.map_view_params = self.view_params_payload(view)
else:
self.geom_data = geom_data
self.property_data = property_data
self.map_view_params = self.view_params_payload(view)
# update current view
cv = self.currentView = copy.deepcopy(view)
# set model ids based on domain
if cv.colorby == 'cell':
self.ids = self.cell_ids
domain = cv.cells
source = self.modelCells
else:
self.ids = self.mat_ids
domain = cv.materials
source = self.modelMaterials
# generate colors if not present
for cell_id, cell in cv.cells.items():
if cell.color is None:
cv.cells.set_color(cell_id, random_rgb())
for mat_id, mat in cv.materials.items():
if mat.color is None:
cv.material.set_color(mat_id, random_rgb())
# construct image data
domain[_OVERLAP] = DomainView(_OVERLAP, "Overlap", cv.overlap_color)
domain[_NOT_FOUND] = DomainView(_NOT_FOUND, "Not Found", cv.domainBackground)
u, inv = np.unique(self.ids, return_inverse=True)
image = np.array([domain[id].color for id in u])[inv]
image.shape = (cv.v_res, cv.h_res, 3)
if cv.masking:
for id, dom in domain.items():
if dom.masked:
image[self.ids == int(id)] = cv.maskBackground
if cv.highlighting:
for id, dom in domain.items():
if dom.highlight:
image[self.ids == int(id)] = cv.highlightBackground
# set model image
self.image = image
# tally data
self.tally_data = None
self.property_data[self.property_data < 0.0] = np.nan
self.temperatures = self.property_data[..., _PROPERTY_INDICES['temperature']]
self.densities = self.property_data[..., _PROPERTY_INDICES['density']]
minmax = {}
for prop in _MODEL_PROPERTIES:
idx = _PROPERTY_INDICES[prop]
prop_data = self.property_data[:, :, idx]
minmax[prop] = (np.min(np.nan_to_num(prop_data)),
np.max(np.nan_to_num(prop_data)))
cv.data_minmax = minmax
if self.activeView.view_params == cv.view_params:
self.activeView.data_minmax = minmax
def undo(self):
""" Revert to previous PlotView instance. Re-generate plot image """
if self.previousViews:
self.subsequentViews.append(copy.deepcopy(self.currentView))
self.activeView = self.previousViews.pop()
def redo(self):
""" Revert to subsequent PlotView instance. Re-generate plot image """
if self.subsequentViews:
self.storeCurrent()
self.activeView = self.subsequentViews.pop()
def getExternalSourceSites(self, n=100):
"""Plot source sites from a source file
"""
if n == 0:
self.source_sites = None
return
sites = openmc.lib.sample_external_source(n)
self.sourceSites = np.array([s.r for s in sites[:n]], dtype=float)
def storeCurrent(self):
""" Add current view to previousViews list """
self.previousViews.append(copy.deepcopy(self.currentView))
def create_tally_image(self, view: Optional[PlotView] = None):
"""
Parameters
----------
view : PlotView
View used to set bounds of the tally data
Returns
-------
tuple
image data (numpy.ndarray), data extents (optional),
data_min_value (float), data_max_value (float),
data label (str)
"""
if view is None:
view = self.currentView
tally_id = view.selectedTally
scores = self.appliedScores
nuclides = self.appliedNuclides
tally_selected = view.selectedTally is not None
tally_visible = view.tallyDataVisible
visible_selection = scores and nuclides
if not tally_selected or not tally_visible or not visible_selection:
return (None, None, None, None, None)
tally = self.statepoint.tallies[tally_id]
tally_value = _TALLY_VALUES[view.tallyValue]
# check score units
units = {_SCORE_UNITS.get(score, _REACTION_UNITS) for score in scores}
if len(units) != 1:
msg_box = QMessageBox()
unit_str = " ".join(units)
msg = "The scores selected have incompatible units:\n"
for unit in units:
msg += " - {}\n".format(unit)
msg_box.setText(msg)
msg_box.setIcon(QMessageBox.Information)
msg_box.setStandardButtons(QMessageBox.Ok)
msg_box.exec()
return (None, None, None, None, None)
units_out = list(units)[0]
contains_distribcell = tally.contains_filter(openmc.DistribcellFilter)
contains_cellinstance = tally.contains_filter(openmc.CellInstanceFilter)
if tally.contains_filter(openmc.MeshFilter):
# Check for volume normalization in order to change units
if view.tallyVolumeNorm:
units_out = _SCORE_UNITS_VOL.get(scores[0], _REACTION_UNITS_VOL)
if tally_value == 'rel_err':
# Get both the std. dev. data and mean data to create the
# relative error data
mean_data = self._create_tally_mesh_image(
tally, 'mean', scores, nuclides, view)[0]
std_dev_data = self._create_tally_mesh_image(
tally, 'std_dev', scores, nuclides, view)[0]
with np.errstate(divide='ignore', invalid='ignore'):
image_data = 100.0 * std_dev_data / mean_data
if np.isnan(image_data).all():
data_min, data_max = 0., 1.
else:
data_min = np.nanmin(image_data)
data_max = np.nanmax(image_data)
return image_data, None, data_min, data_max, '% error'
else:
image = self._create_tally_mesh_image(tally,
tally_value,
scores,
nuclides,
view)
return image + (units_out,)
elif tally.contains_filter(openmc.MeshMaterialFilter):
image = self._create_tally_filter_image(
tally, tally_value, openmc.MeshMaterialFilter, scores, nuclides, view)
return image + (units_out,)
elif contains_distribcell or contains_cellinstance:
if tally_value == 'rel_err':
mean_data = self._create_distribcell_image(
tally, 'mean', scores, nuclides, contains_cellinstance)
std_dev_data = self._create_distribcell_image(
tally, 'std_dev', scores, nuclides)
image_data = 100 * np.divide(
std_dev_data[0], mean_data[0],
out=np.zeros_like(mean_data[0]),
where=mean_data != 0
)
data_min = np.min(image_data)
data_max = np.max(image_data)
return image_data, None, data_min, data_max, '% error'
else:
image = self._create_distribcell_image(
tally, tally_value, scores, nuclides, contains_cellinstance)
return image + (units_out,)
else:
# same as above, get the std. dev. data
# and mean date to produce the relative error data
if tally_value == 'rel_err':
mean_data = self._create_tally_domain_image(tally,
'mean',
scores,
nuclides)
std_dev_data = self._create_tally_domain_image(tally,
'std_dev',
scores,
nuclides)
image_data = 100 * np.divide(std_dev_data[0],
mean_data[0],
out=np.zeros_like(mean_data[0]),
where=mean_data != 0)
# adjust for NaNs in bins without tallies
image_data = np.nan_to_num(image_data,
nan=0.0,
posinf=0.0,
neginf=0.0)
extents = mean_data[1]
data_min = np.min(image_data)
data_max = np.max(image_data)
return image_data, extents, data_min, data_max, '% error'
else:
image = self._create_tally_domain_image(tally,
tally_value,
scores,
nuclides)
return image + (units_out,)
def _create_tally_domain_image(self, tally, tally_value, scores, nuclides, view=None):
# data resources used throughout
if view is None:
view = self.currentView
data = tally.get_reshaped_data(tally_value)
data_out = np.full(self.ids.shape, -1.0)
def _do_op(array, tally_value, ax=0):
if tally_value == 'mean':
return np.sum(array, axis=ax)
elif tally_value == 'std_dev':
return np.sqrt(np.sum(array**2, axis=ax))
# data structure for tracking which spatial
# filter bins are enabled
spatial_filter_bins = defaultdict(list)
n_spatial_filters = 0
for tally_filter in tally.filters:
if tally_filter in self.appliedFilters:
selected_bins = self.appliedFilters[tally_filter]
if type(tally_filter) in _SPATIAL_FILTERS:
spatial_filter_bins[tally_filter] = selected_bins
n_spatial_filters += 1
else:
slc = [slice(None)] * len(data.shape)
slc[n_spatial_filters] = selected_bins
slc = tuple(slc)
data = _do_op(data[slc], tally_value, n_spatial_filters)
else:
data[:] = 0.0
data = _do_op(data, tally_value, n_spatial_filters)
# filter by selected scores
selected_scores = []
for idx, score in enumerate(tally.scores):
if score in scores:
selected_scores.append(idx)
data = _do_op(data[..., np.array(selected_scores)], tally_value, -1)
# filter by selected nuclides
selected_nuclides = []
for idx, nuclide in enumerate(tally.nuclides):
if nuclide in nuclides:
selected_nuclides.append(idx)
data = _do_op(data[..., np.array(selected_nuclides)], tally_value, -1)
# get data limits
data_min = np.min(data)
data_max = np.max(data)
# for all combinations of spatial bins, create a mask
# and set image data values
spatial_filters = list(spatial_filter_bins.keys())
spatial_bins = list(spatial_filter_bins.values())
for bin_indices in itertools.product(*spatial_bins):
# look up the tally value
tally_val = data[bin_indices]
if tally_val == 0.0:
continue
# generate a mask with the correct size
mask = np.full(self.ids.shape, True, dtype=bool)
for tally_filter, bin_idx in zip(spatial_filters, bin_indices):
bin = tally_filter.bins[bin_idx]
if isinstance(tally_filter, openmc.CellFilter):
mask &= self.cell_ids == bin
elif isinstance(tally_filter, openmc.MaterialFilter):
mask &= self.mat_ids == bin
elif isinstance(tally_filter, openmc.UniverseFilter):
# get the statepoint summary
univ_cells = self.statepoint.universes[bin].cells
for cell in univ_cells:
mask &= self.cell_ids == cell
# set image data values
data_out[mask] = tally_val
# mask out invalid values
image_data = np.ma.masked_where(data_out < 0.0, data_out)
return image_data, None, data_min, data_max
def _create_distribcell_image(self, tally, tally_value, scores, nuclides, cellinstance=False):
# Get flattened array of tally results
data = tally.get_values(scores=scores, nuclides=nuclides, value=tally_value)
data = data.flatten()
# Create an empty array of appropriate shape for image
image_data = np.full_like(self.ids, np.nan, dtype=float)
# Determine mapping of cell IDs to list of (instance, tally value).
if cellinstance:
f = tally.find_filter(openmc.CellInstanceFilter)
cell_id_to_inst_value = defaultdict(list)
for value, (cell_id, instance) in zip(data, f.bins):
cell_id_to_inst_value[cell_id].append((instance, value))
else:
f = tally.find_filter(openmc.DistribcellFilter)
cell_id_to_inst_value = {f.bins[0]: list(enumerate(data))}
for cell_id, value_list in cell_id_to_inst_value.items():
# Get mask for each relevant cell
cell_id_mask = (self.cell_ids == cell_id)
# For each cell, iterate over instances and corresponding tally
# values and set any matching pixels
for instance, value in value_list:
instance_mask = (self.instances == instance)
image_data[cell_id_mask & instance_mask] = value
data_min = np.min(data)
data_max = np.max(data)
image_data = np.ma.masked_where(image_data < 0.0, image_data)
return image_data, None, data_min, data_max
def cpp_mesh_ids(self):
return list(openmc.lib.meshes.keys())
def mesh_plot_bins(self, mesh_id, view: PlotView = None, translation: tuple[float, float, float] = None):
mesh = openmc.lib.meshes[mesh_id]
if view is None:
view = self.currentView
if translation is None:
origin = view.origin
else:
origin = (view.origin[0] - translation[0],
view.origin[1] - translation[1],
view.origin[2] - translation[2])
mesh_bins = mesh.get_plot_bins(
origin=origin,
width=(view.width, view.height),
basis=view.basis,
pixels=(view.h_res, view.v_res),
)
return mesh_bins
def _create_tally_mesh_image(
self, tally: openmc.Tally, tally_value: TallyValueType,
scores: Tuple[str], nuclides: Tuple[str], view: PlotView = None
):
# some variables used throughout
if view is None:
view = self.currentView
mesh_filter = tally.find_filter(openmc.MeshFilter)
mesh = mesh_filter.mesh
def _do_op(array, tally_value, ax=0):
if tally_value == 'mean':
return np.sum(array, axis=ax)
elif tally_value == 'std_dev':
return np.sqrt(np.sum(array**2, axis=ax))
# start with reshaped data
data = tally.get_reshaped_data(tally_value)
# move mesh axes to the end of the filters
filter_idx = [type(filter) for filter in tally.filters].index(openmc.MeshFilter)
data = np.moveaxis(data, filter_idx, -1)
# sum over the rest of the tally filters
for tally_filter in tally.filters:
if type(tally_filter) == openmc.MeshFilter:
continue
selected_bins = self.appliedFilters[tally_filter]
if selected_bins:
# sum filter data for the selected bins
data = data[np.array(selected_bins)].sum(axis=0)
else:
# if the filter is completely unselected,
# set all of its data to zero and remove the axis
data[:] = 0.0
data = _do_op(data, tally_value)
# filter by selected nuclides
if not nuclides:
data = 0.0
selected_nuclides = []
for idx, nuclide in enumerate(tally.nuclides):
if nuclide in nuclides:
selected_nuclides.append(idx)
data = _do_op(data[np.array(selected_nuclides)], tally_value)
# filter by selected scores
if not scores:
data = 0.0
selected_scores = []
for idx, score in enumerate(tally.scores):
if score in scores:
selected_scores.append(idx)
data = _do_op(data[np.array(selected_scores)], tally_value)
# Get mesh bins from openmc.lib
mesh_cpp = openmc.lib.meshes[mesh.id]
mesh_bins = self.mesh_plot_bins(mesh.id, view, mesh_filter.translation)
# Apply volume normalization
if view.tallyVolumeNorm:
data /= mesh_cpp.volumes
# set image data
image_data = np.full_like(self.ids, np.nan, dtype=float)
mask = (mesh_bins >= 0)
image_data[mask] = data[mesh_bins[mask]]
# get dataset's min/max
data_min = np.min(data)
data_max = np.max(data)
return image_data, None, data_min, data_max
def _create_tally_filter_image(
self, tally: openmc.Tally, tally_value: TallyValueType, filter_class,
scores: Tuple[str], nuclides: Tuple[str], view: PlotView = None
):
# some variables used throughout
if view is None:
view = self.currentView
def _do_op(array, tally_value, ax=0):
if tally_value == 'mean':
return np.sum(array, axis=ax)
elif tally_value == 'std_dev':
return np.sqrt(np.sum(array**2, axis=ax))
# start with reshaped data