openmc/openmc/plotter.py at develop · openmc-dev/openmc · GitHub

History

985 lines (872 loc) · 40.3 KB

Raw

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

from __future__ import annotations

from itertools import chain

from numbers import Integral, Real

from typing import Dict, Iterable, List

import numpy as np

import openmc.checkvalue as cv

import openmc.data

# Supported keywords for continuous-energy cross section plotting

PLOT_TYPES = {'total', 'scatter', 'elastic', 'inelastic', 'fission',

'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity',

'slowing-down power', 'damage'}

# Supported keywords for multi-group cross section plotting

PLOT_TYPES_MGXS = {'total', 'absorption', 'scatter', 'fission',

'kappa-fission', 'nu-fission', 'prompt-nu-fission',

'deleyed-nu-fission', 'chi', 'chi-prompt', 'chi-delayed',

'inverse-velocity', 'beta', 'decay-rate', 'unity'}

# Create a dictionary which can be used to convert PLOT_TYPES_MGXS to the

# openmc.XSdata attribute name needed to access the data

_PLOT_MGXS_ATTR = {line: line.replace(' ', '_').replace('-', '_')

for line in PLOT_TYPES_MGXS}

_PLOT_MGXS_ATTR['scatter'] = 'scatter_matrix'

# Special MT values

UNITY_MT = -1

XI_MT = -2

# MTs to combine to generate associated plot_types

_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt != 27]

PLOT_TYPES_MT = {

'total': openmc.data.SUM_RULES[1],

'scatter': [2] + _INELASTIC,

'elastic': [2],

'inelastic': _INELASTIC,

'fission': [18],

'absorption': [27],

'capture': [101],

'nu-fission': [18],

'nu-scatter': [2] + _INELASTIC,

'unity': [UNITY_MT],

'slowing-down power': [2] + [XI_MT],

'damage': [444]

}

# Types of plots to plot linearly in y

PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter',

'nu-fission / absorption', 'fission / absorption'}

# Minimum and maximum energies for plotting (units of eV)

_MIN_E = 1.e-5

_MAX_E = 20.e6

ELEMENT_NAMES = list(openmc.data.ELEMENT_SYMBOL.values())[1:]

def _get_legend_label(this, type):

"""Gets a label for the element or nuclide or material and reaction plotted"""

if isinstance(this, str):

if type in openmc.data.DADZ:

if this in ELEMENT_NAMES:

return f'{this} {type}'

else: # this is a nuclide so the legend can contain more information

z, a, m = openmc.data.zam(this)

da, dz = openmc.data.DADZ[type]

gnds_name = openmc.data.gnds_name(z + dz, a + da, m)

# makes a string with nuclide reaction and new nuclide

# For example "Be9 (n,2n) Be8"

return f'{this} {type} {gnds_name}'

return f'{this} {type}'

elif this.name == '':

return f'Material {this.id} {type}'

else:

return f'{this.name} {type}'

def _get_yaxis_label(reactions, divisor_types):

"""Gets a y axis label for the type of data plotted"""

heat_values = {"heating", "heating-local", "damage-energy"}

# if all the types are heating a different stem and unit is needed

if all(set(value).issubset(heat_values) for value in reactions.values()):

stem = "Heating"

elif all(isinstance(item, str) for item in reactions.keys()):

for nuc_reactions in reactions.values():

for reaction in nuc_reactions:

if reaction in heat_values:

raise TypeError(

"Mixture of heating and Microscopic reactions. "

"Invalid type for plotting"

)

stem = "Microscopic"

elif all(isinstance(item, openmc.Material) for item in reactions.keys()):

stem = 'Macroscopic'

else:

msg = "Mixture of openmc.Material and elements/nuclides. Invalid type for plotting"

raise TypeError(msg)

if divisor_types:

mid, units = "Data", ""

else:

mid = "Cross Section"

units = {

"Macroscopic": "[1/cm]",

"Microscopic": "[b]",

"Heating": "[eV-barn]",

}[stem]

return f'{stem} {mid} {units}'

def _get_title(reactions):

"""Gets a title for the type of data plotted"""

if len(reactions) == 1:

this, = reactions

name = this.name if isinstance(this, openmc.Material) else this

return f'Cross Section Plot For {name}'

else:

return 'Cross Section Plot'

def plot_xs(

reactions: Dict[str | openmc.Material, List[str]],

divisor_types: Iterable[str] | None = None,

temperature: float = 294.0,

axis: "plt.Axes" | None = None,

sab_name: str | None = None,

ce_cross_sections: str | None = None,

mg_cross_sections: str | None = None,

enrichment: float | None = None,

plot_CE: bool = True,

orders: Iterable[int] | None = None,

divisor_orders: Iterable[int] | None = None,

energy_axis_units: str = "eV",

**kwargs,

) -> "plt.Figure" | None:

"""Creates a figure of continuous-energy cross sections for this item.

Parameters

----------

reactions : dict

keys can be either a nuclide or element in string form or an

openmc.Material object. Values are a list of the types of

cross sections to include in the plot.

divisor_types : Iterable of values of PLOT_TYPES, optional

Cross section types which will divide those produced by types

before plotting. A type of 'unity' can be used to effectively not

divide some types.

temperature : float, optional

Temperature in Kelvin to plot. If not specified, a default

temperature of 294K will be plotted. Note that the nearest

temperature in the library for each nuclide will be used as opposed

to using any interpolation.

axis : matplotlib.axes, optional

A previously generated axis to use for plotting. If not specified,

a new axis and figure will be generated.

sab_name : str, optional

Name of S(a,b) library to apply to MT=2 data when applicable.

ce_cross_sections : str, optional

Location of cross_sections.xml file. Default is None.

mg_cross_sections : str, optional

Location of MGXS HDF5 Library file. Default is None.

enrichment : float, optional

Enrichment for U235 in weight percent. For example, input 4.95 for

4.95 weight percent enriched U. Default is None.

plot_CE : bool, optional

Denotes whether or not continuous-energy will be plotted. Defaults to

plotting the continuous-energy data.

orders : Iterable of Integral, optional

The scattering order or delayed group index to use for the

corresponding entry in types. Defaults to the 0th order for scattering

and the total delayed neutron data. This only applies to plots of

multi-group data.

divisor_orders : Iterable of Integral, optional

Same as orders, but for divisor_types

**kwargs :

All keyword arguments are passed to

:func:`matplotlib.pyplot.figure`.

energy_axis_units : {'eV', 'keV', 'MeV'}

Units used on the plot energy axis

.. versionadded:: 0.15.0

Returns

-------

fig : matplotlib.figure.Figure

If axis is None, then a Matplotlib Figure of the generated

cross section will be returned. Otherwise, a value of

None will be returned as the figure and axes have already been

generated.

"""

import matplotlib.pyplot as plt

cv.check_type("plot_CE", plot_CE, bool)

cv.check_value("energy_axis_units", energy_axis_units, {"eV", "keV", "MeV"})

axis_scaling_factor = {"eV": 1.0, "keV": 1e-3, "MeV": 1e-6}

# Generate the plot

if axis is None:

fig, ax = plt.subplots(**kwargs)

else:

fig = None

ax = axis

all_types = []

for this, types in reactions.items():

all_types = all_types + types

if plot_CE:

cv.check_type("this", this, (str, openmc.Material))

# Calculate for the CE cross sections

E, data = calculate_cexs(this, types, temperature, sab_name,

ce_cross_sections, enrichment)

if divisor_types:

cv.check_length('divisor types', divisor_types, len(types))

Ediv, data_div = calculate_cexs(this, divisor_types, temperature,

sab_name, ce_cross_sections,

enrichment)

# Create a new union grid, interpolate data and data_div on to that

# grid, and then do the actual division

Enum = E[:]

E = np.union1d(Enum, Ediv)

data_new = np.zeros((len(types), len(E)))

for line in range(len(types)):

data_new[line, :] = \

np.divide(np.interp(E, Enum, data[line, :]),

np.interp(E, Ediv, data_div[line, :]))

if divisor_types[line] != 'unity':

types[line] = types[line] + ' / ' + divisor_types[line]

data = data_new

else:

# Calculate for MG cross sections

E, data = calculate_mgxs(this, types, orders, temperature,

mg_cross_sections, ce_cross_sections,

enrichment)

if divisor_types:

cv.check_length('divisor types', divisor_types, len(types))

Ediv, data_div = calculate_mgxs(this, divisor_types,

divisor_orders, temperature,

mg_cross_sections,

ce_cross_sections, enrichment)

# Perform the division

for line in range(len(types)):

data[line, :] /= data_div[line, :]

if divisor_types[line] != 'unity':

types[line] += ' / ' + divisor_types[line]

E *= axis_scaling_factor[energy_axis_units]

# Plot the data

for i in range(len(data)):

data[i, :] = np.nan_to_num(data[i, :])

if np.sum(data[i, :]) > 0.:

ax.plot(E, data[i, :], label=_get_legend_label(this, types[i]))

# Set to loglog or semilogx depending on if we are plotting a data

# type which we expect to vary linearly

if set(all_types).issubset(PLOT_TYPES_LINEAR):

ax.set_xscale('log')

ax.set_yscale('linear')

else:

ax.set_xscale('log')

ax.set_yscale('log')

ax.set_xlabel(f"Energy [{energy_axis_units}]")

if plot_CE:

ax.set_xlim(

_MIN_E * axis_scaling_factor[energy_axis_units],

_MAX_E * axis_scaling_factor[energy_axis_units],

)

else:

ax.set_xlim(E[-1], E[0])

ax.set_ylabel(_get_yaxis_label(reactions, divisor_types))

ax.legend(loc='best')

ax.set_title(_get_title(reactions))

return fig

def calculate_cexs(this, types, temperature=294., sab_name=None,

cross_sections=None, enrichment=None, ncrystal_cfg=None):

"""Calculates continuous-energy cross sections of a requested type.

Parameters

----------

this : str or openmc.Material

Object to source data from. Nuclides and elements should be input as a

str

types : Iterable of values of PLOT_TYPES

The type of cross sections to calculate

temperature : float, optional

Temperature in Kelvin to plot. If not specified, a default

temperature of 294K will be plotted. Note that the nearest

temperature in the library for each nuclide will be used as opposed

to using any interpolation.

sab_name : str, optional

Name of S(a,b) library to apply to MT=2 data when applicable.

cross_sections : str, optional

Location of cross_sections.xml file. Default is None.

enrichment : float, optional

Enrichment for U235 in weight percent. For example, input 4.95 for

4.95 weight percent enriched U. Default is None

(natural composition).

ncrystal_cfg : str, optional

Configuration string for NCrystal material.

Returns

-------

energy_grid : numpy.ndarray

Energies at which cross sections are calculated, in units of eV

data : numpy.ndarray

Cross sections calculated at the energy grid described by energy_grid

"""

# Check types

cv.check_type('this', this, (str, openmc.Material))

cv.check_type('temperature', temperature, Real)

if sab_name:

cv.check_type('sab_name', sab_name, str)

if enrichment:

cv.check_type('enrichment', enrichment, Real)

if isinstance(this, str):

if this in ELEMENT_NAMES:

energy_grid, data = _calculate_cexs_elem_mat(

this, types, temperature, cross_sections, sab_name, enrichment

)

else:

energy_grid, xs = _calculate_cexs_nuclide(

this, types, temperature, sab_name, cross_sections,

ncrystal_cfg

)

# Convert xs (Iterable of Callable) to a grid of cross section values

# calculated on the points in energy_grid for consistency with the

# element and material functions.

data = np.zeros((len(types), len(energy_grid)))

for line in range(len(types)):

data[line, :] = xs[line](energy_grid)

else:

energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature,

cross_sections)

return energy_grid, data

def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None,

cross_sections=None, ncrystal_cfg=None):

"""Calculates continuous-energy cross sections of a requested type.

Parameters

----------

this : str

Nuclide object to source data from

types : Iterable of str or Integral

The type of cross sections to calculate; values can either be those

in openmc.PLOT_TYPES or keys from openmc.data.REACTION_MT which

correspond to a reaction description e.g '(n,2n)' or integers which

correspond to reaction channel (MT) numbers.

temperature : float, optional

Temperature in Kelvin to plot. If not specified, a default

temperature of 294K will be plotted. Note that the nearest

temperature in the library for each nuclide will be used as opposed

to using any interpolation.

sab_name : str, optional

Name of S(a,b) library to apply to MT=2 data when applicable.

cross_sections : str, optional

Location of cross_sections.xml file. Default is None.

ncrystal_cfg : str, optional

Configuration string for NCrystal material.

Returns

-------

energy_grid : numpy.ndarray

Energies at which cross sections are calculated, in units of eV

data : Iterable of Callable

Requested cross section functions

"""

# Load the library

library = openmc.data.DataLibrary.from_xml(cross_sections)

# Convert temperature to format needed for access in the library

strT = f"{int(round(temperature))}K"

T = temperature

# Now we can create the data sets to be plotted

energy_grid = []

xs = []

lib = library.get_by_material(this)

if lib is not None:

nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path'])

# Obtain the nearest temperature

if strT in nuc.temperatures:

nucT = strT

else:

delta_T = np.array(nuc.kTs) - T * openmc.data.K_BOLTZMANN

closest_index = np.argmin(np.abs(delta_T))

nucT = nuc.temperatures[closest_index]

# Prep S(a,b) data if needed

if sab_name:

sab = openmc.data.ThermalScattering.from_hdf5(

library.get_by_material(sab_name, data_type='thermal')['path'])

# Obtain the nearest temperature

if strT in sab.temperatures:

sabT = strT

else:

delta_T = np.array(sab.kTs) - T * openmc.data.K_BOLTZMANN

closest_index = np.argmin(np.abs(delta_T))

sabT = sab.temperatures[closest_index]

# Create an energy grid composed the S(a,b) and the nuclide's grid

grid = nuc.energy[nucT]

sab_Emax = 0.

sab_funcs = []

if sab.elastic is not None:

elastic = sab.elastic.xs[sabT]

if isinstance(elastic, openmc.data.CoherentElastic):

grid = np.union1d(grid, elastic.bragg_edges)

if elastic.bragg_edges[-1] > sab_Emax:

sab_Emax = elastic.bragg_edges[-1]

elif isinstance(elastic, openmc.data.Tabulated1D):

grid = np.union1d(grid, elastic.x)

if elastic.x[-1] > sab_Emax:

sab_Emax = elastic.x[-1]

sab_funcs.append(elastic)

if sab.inelastic is not None:

inelastic = sab.inelastic.xs[sabT]

grid = np.union1d(grid, inelastic.x)

if inelastic.x[-1] > sab_Emax:

sab_Emax = inelastic.x[-1]

sab_funcs.append(inelastic)

energy_grid = grid

else:

energy_grid = nuc.energy[nucT]

# Parse the types

mts = []

ops = []

yields = []

for line in types:

if line in PLOT_TYPES:

tmp_mts = [mtj for mti in PLOT_TYPES_MT[line] for mtj in

nuc.get_reaction_components(mti)]

mts.append(tmp_mts)

if line.startswith('nu'):

yields.append(True)

else:

yields.append(False)

if XI_MT in tmp_mts:

ops.append((np.add,) * (len(tmp_mts) - 2) + (np.multiply,))

else:

ops.append((np.add,) * (len(tmp_mts) - 1))

elif line in openmc.data.REACTION_MT:

mt_number = openmc.data.REACTION_MT[line]

cv.check_type('MT in types', mt_number, Integral)

cv.check_greater_than('MT in types', mt_number, 0)

tmp_mts = nuc.get_reaction_components(mt_number)

mts.append(tmp_mts)

ops.append((np.add,) * (len(tmp_mts) - 1))

yields.append(False)

elif isinstance(line, int):

# Not a built-in type, we have to parse it ourselves

cv.check_type('MT in types', line, Integral)

cv.check_greater_than('MT in types', line, 0)

tmp_mts = nuc.get_reaction_components(line)

mts.append(tmp_mts)

ops.append((np.add,) * (len(tmp_mts) - 1))

yields.append(False)

else:

raise TypeError("Invalid type", line)

for i, mt_set in enumerate(mts):

# Get the reaction xs data from the nuclide

funcs = []

op = ops[i]

for mt in mt_set:

if mt == 2:

if sab_name:

# Then we need to do a piece-wise function of

# The S(a,b) and non-thermal data

sab_sum = openmc.data.Sum(sab_funcs)

pw_funcs = openmc.data.Regions1D(

[sab_sum, nuc[mt].xs[nucT]],

[sab_Emax])

funcs.append(pw_funcs)

elif ncrystal_cfg:

import NCrystal

nc_scatter = NCrystal.createScatter(ncrystal_cfg)

nc_func = nc_scatter.xsect

nc_emax = 5 # eV # this should be obtained from NCRYSTAL_MAX_ENERGY

energy_grid = np.union1d(np.geomspace(min(energy_grid),

1.1*nc_emax,

1000),energy_grid) # NCrystal does not have

# an intrinsic energy grid

pw_funcs = openmc.data.Regions1D(

[nc_func, nuc[mt].xs[nucT]],

[nc_emax])

funcs.append(pw_funcs)

else:

funcs.append(nuc[mt].xs[nucT])

elif mt in nuc:

if yields[i]:

# Get the total yield first if available. This will be

# used primarily for fission.

for prod in chain(nuc[mt].products,

nuc[mt].derived_products):

if prod.particle == 'neutron' and \

prod.emission_mode == 'total':

func = openmc.data.Combination(

[nuc[mt].xs[nucT], prod.yield_],

[np.multiply])

funcs.append(func)

break

else:

# Total doesn't exist so we have to create from

# prompt and delayed. This is used for scatter

# multiplication.

func = None

for prod in chain(nuc[mt].products,

nuc[mt].derived_products):

if prod.particle == 'neutron' and \

prod.emission_mode != 'total':

if func:

func = openmc.data.Combination(

[prod.yield_, func], [np.add])

else:

func = prod.yield_

if func:

funcs.append(openmc.data.Combination(

[func, nuc[mt].xs[nucT]], [np.multiply]))

else:

# If func is still None, then there were no

# products. In that case, assume the yield is

# one as its not provided for some summed

# reactions like MT=4

funcs.append(nuc[mt].xs[nucT])

else:

funcs.append(nuc[mt].xs[nucT])

elif mt == UNITY_MT:

funcs.append(lambda x: 1.)

elif mt == XI_MT:

awr = nuc.atomic_weight_ratio

alpha = ((awr - 1.) / (awr + 1.))**2

xi = 1. + alpha * np.log(alpha) / (1. - alpha)

funcs.append(lambda x: xi)

else:

funcs.append(lambda x: 0.)

funcs = funcs if funcs else [lambda x: 0.]

xs.append(openmc.data.Combination(funcs, op))

else:

raise ValueError(this + " not in library")

return energy_grid, xs

def _calculate_cexs_elem_mat(this, types, temperature=294.,

cross_sections=None, sab_name=None,

enrichment=None):

"""Calculates continuous-energy cross sections of a requested type.

Parameters

----------

this : openmc.Material or str

Object to source data from. Element can be input as str

types : Iterable of values of PLOT_TYPES

The type of cross sections to calculate

temperature : float, optional

Temperature in Kelvin to plot. If not specified, a default

temperature of 294K will be plotted. Note that the nearest

temperature in the library for each nuclide will be used as opposed

to using any interpolation.

cross_sections : str, optional

Location of cross_sections.xml file. Default is None.

sab_name : str, optional

Name of S(a,b) library to apply to MT=2 data when applicable.

enrichment : float, optional

Enrichment for U235 in weight percent. For example, input 4.95 for

4.95 weight percent enriched U. Default is None

(natural composition).

Returns

-------

energy_grid : numpy.ndarray

Energies at which cross sections are calculated, in units of eV

data : numpy.ndarray

Cross sections calculated at the energy grid described by energy_grid

"""

if isinstance(this, openmc.Material):

if this.temperature is not None:

T = this.temperature

else:

T = temperature

else:

T = temperature

# Load the library

library = openmc.data.DataLibrary.from_xml(cross_sections)

ncrystal_cfg = None

if isinstance(this, openmc.Material):

# Expand elements in to nuclides with atomic densities

nuc_fractions = this.get_nuclide_atom_densities()

# Create a dict of [nuclide name] = nuclide object to carry forward

# with a common nuclides format between openmc.Material and Elements

nuclides = {nuclide: nuclide for nuclide in nuc_fractions}

# Add NCrystal cfg string if it exists

ncrystal_cfg = this.ncrystal_cfg

else:

# Expand elements in to nuclides with atomic densities

nuclides = openmc.Element(this).expand(1., 'ao', enrichment=enrichment,

cross_sections=cross_sections)

# For ease of processing split out the nuclide and its fraction

nuc_fractions = {nuclide[0]: nuclide[1] for nuclide in nuclides}

# Create a dict of [nuclide name] = nuclide object to carry forward

# with a common nuclides format between openmc.Material and Elements

nuclides = {nuclide[0]: nuclide[0] for nuclide in nuclides}

# Identify the nuclides which have S(a,b) data

sabs = {}

for nuclide in nuclides.items():

sabs[nuclide[0]] = None

if isinstance(this, openmc.Material):

for sab_name, _ in this._sab:

sab = openmc.data.ThermalScattering.from_hdf5(

library.get_by_material(sab_name, data_type='thermal')['path'])

for nuc in sab.nuclides:

sabs[nuc] = sab_name

else:

if sab_name:

sab = openmc.data.ThermalScattering.from_hdf5(

library.get_by_material(sab_name, data_type='thermal')['path'])

for nuc in sab.nuclides:

sabs[nuc] = sab_name

# Now we can create the data sets to be plotted

xs = {}

E = []

for nuclide in nuclides.items():

name = nuclide[0]

nuc = nuclide[1]

sab_name = sabs[name]

temp_E, temp_xs = calculate_cexs(nuc, types, T, sab_name, cross_sections,

ncrystal_cfg=ncrystal_cfg

)

E.append(temp_E)

# Since the energy grids are different, store the cross sections as

# a tabulated function so they can be calculated on any grid needed.

xs[name] = [openmc.data.Tabulated1D(temp_E, temp_xs[line])

for line in range(len(types))]

# Condense the data for every nuclide

# First create a union energy grid

energy_grid = E[0]

for grid in E[1:]:

energy_grid = np.union1d(energy_grid, grid)

# Now we can combine all the nuclidic data

data = np.zeros((len(types), len(energy_grid)))

for line in range(len(types)):

if types[line] == 'unity':

data[line, :] = 1.

else:

for nuclide in nuclides.items():

name = nuclide[0]

data[line, :] += (nuc_fractions[name] *

xs[name][line](energy_grid))

return energy_grid, data

def calculate_mgxs(this, types, orders=None, temperature=294.,

cross_sections=None, ce_cross_sections=None,

enrichment=None):

"""Calculates multi-group cross sections of a requested type.

If the data for the nuclide or macroscopic object in the library is

represented as angle-dependent data then this method will return the

geometric average cross section over all angles.

Parameters

----------

this : str or openmc.Material

Object to source data from. Nuclides and elements can be input as a str

types : Iterable of values of PLOT_TYPES_MGXS

The type of cross sections to calculate

orders : Iterable of Integral, optional

The scattering order or delayed group index to use for the

corresponding entry in types. Defaults to the 0th order for scattering

and the total delayed neutron data.

temperature : float, optional

Temperature in Kelvin to plot. If not specified, a default

temperature of 294K will be plotted. Note that the nearest

temperature in the library for each nuclide will be used as opposed

to using any interpolation.

cross_sections : str, optional

Location of MGXS HDF5 Library file. Default is None.

ce_cross_sections : str, optional

Location of continuous-energy cross_sections.xml file. Default is None.

enrichment : float, optional

Enrichment for U235 in weight percent. For example, input 4.95 for

4.95 weight percent enriched U. Default is None

(natural composition).

Returns

-------

energy_grid : numpy.ndarray

Energies at which cross sections are calculated, in units of eV

data : numpy.ndarray

Cross sections calculated at the energy grid described by energy_grid

"""

# Check types

cv.check_type('temperature', temperature, Real)

if enrichment:

cv.check_type('enrichment', enrichment, Real)

cv.check_iterable_type('types', types, str)

cv.check_type("cross_sections", cross_sections, str)

library = openmc.MGXSLibrary.from_hdf5(cross_sections)

if this in ELEMENT_NAMES or isinstance(this, openmc.Material):

mgxs = _calculate_mgxs_elem_mat(this, types, library, orders,

temperature, ce_cross_sections,

enrichment)

elif isinstance(this, str):

mgxs = _calculate_mgxs_nuc_macro(this, types, library, orders,

temperature)

else:

raise TypeError("Invalid type")

# Convert the data to the format needed

data = np.zeros((len(types), 2 * library.energy_groups.num_groups))

energy_grid = np.zeros(2 * library.energy_groups.num_groups)

for g in range(library.energy_groups.num_groups):

energy_grid[g * 2: g * 2 + 2] = \

library.energy_groups.group_edges[g: g + 2]

# Ensure the energy will show on a log-axis by replacing 0s with a

# sufficiently small number

energy_grid[0] = max(energy_grid[0], _MIN_E)

for line in range(len(types)):

for g in range(library.energy_groups.num_groups):

data[line, g * 2: g * 2 + 2] = mgxs[line, g]

return energy_grid[::-1], data

def _calculate_mgxs_nuc_macro(this, types, library, orders=None,

temperature=294.):

"""Determines the multi-group cross sections of a nuclide or macroscopic

object.

If the data for the nuclide or macroscopic object in the library is

represented as angle-dependent data then this method will return the

geometric average cross section over all angles.

Parameters

----------

this : str

Object to source data from

types : Iterable of str

The type of cross sections to calculate; values can either be those

in openmc.PLOT_TYPES_MGXS

library : openmc.MGXSLibrary

MGXS Library containing the data of interest

orders : Iterable of Integral, optional

The scattering order or delayed group index to use for the

corresponding entry in types. Defaults to the 0th order for scattering

and the total delayed neutron data.

temperature : float, optional

Temperature in Kelvin to plot. If not specified, a default

temperature of 294K will be plotted. Note that the nearest

temperature in the library for each nuclide will be used as opposed

to using any interpolation.

Returns

-------

data : numpy.ndarray

Cross sections calculated at the energy grid described by energy_grid

"""

# Check the parameters and grab order/delayed groups

if orders:

cv.check_iterable_type('orders', orders, Integral,

min_depth=len(types), max_depth=len(types))

else:

orders = [None] * len(types)

for i, line in enumerate(types):

cv.check_type("line", line, str)

cv.check_value("line", line, PLOT_TYPES_MGXS)

if orders[i]:

cv.check_greater_than("order value", orders[i], 0, equality=True)

xsdata = library.get_by_name(this)

if xsdata is not None:

# Obtain the nearest temperature

t = np.abs(xsdata.temperatures - temperature).argmin()

# Get the data

data = np.zeros((len(types), library.energy_groups.num_groups))

for i, line in enumerate(types):

if 'fission' in line and not xsdata.fissionable:

continue

elif line == 'unity':

data[i, :] = 1.

else:

# Now we have to get the cross section data and properly

# treat it depending on the requested type.

# First get the data in a generic fashion

temp_data = getattr(xsdata, _PLOT_MGXS_ATTR[line])[t]

shape = temp_data.shape[:]

# If we have angular data, then want the geometric

# average over all provided angles. Since the angles are

# equi-distant, un-weighted averaging will suffice

if xsdata.representation == 'angle':

temp_data = np.mean(temp_data, axis=(0, 1))

# Now we can look at the shape of the data to identify how

# it should be modified to produce an array of values

# with groups.

if shape in (xsdata.xs_shapes["[G']"],

xsdata.xs_shapes["[G]"]):

# Then the data is already an array vs groups so copy

# and move along

data[i, :] = temp_data

elif shape == xsdata.xs_shapes["[G][G']"]:

# Sum the data over outgoing groups to create our array vs

# groups

data[i, :] = np.sum(temp_data, axis=1)

elif shape == xsdata.xs_shapes["[DG]"]:

# Then we have a constant vs groups with a value for each

# delayed group. The user-provided value of orders tells us

# which delayed group we want. If none are provided, then

# we sum all the delayed groups together.

if orders[i]:

if orders[i] < len(shape[0]):

data[i, :] = temp_data[orders[i]]

else:

data[i, :] = np.sum(temp_data[:])

elif shape in (xsdata.xs_shapes["[DG][G']"],

xsdata.xs_shapes["[DG][G]"]):

# Then we have an array vs groups with values for each

# delayed group. The user-provided value of orders tells us

# which delayed group we want. If none are provided, then

# we sum all the delayed groups together.

if orders[i]:

if orders[i] < len(shape[0]):

data[i, :] = temp_data[orders[i], :]

else:

data[i, :] = np.sum(temp_data[:, :], axis=0)

elif shape == xsdata.xs_shapes["[DG][G][G']"]:

# Then we have a delayed group matrix. We will first

# remove the outgoing group dependency

temp_data = np.sum(temp_data, axis=-1)

# And then proceed in exactly the same manner as the

# "[DG][G']" or "[DG][G]" shapes in the previous block.

if orders[i]:

if orders[i] < len(shape[0]):

data[i, :] = temp_data[orders[i], :]

else:

data[i, :] = np.sum(temp_data[:, :], axis=0)

elif shape == xsdata.xs_shapes["[G][G'][Order]"]:

# This is a scattering matrix with angular data

# First remove the outgoing group dependence

temp_data = np.sum(temp_data, axis=1)

# The user either provided a specific order or we resort

# to the default 0th order

if orders[i]:

order = orders[i]

else:

order = 0

# If the order is available, store the data for that order

# if it is not available, then the expansion coefficient

# is zero and thus we already have the correct value.

if order < shape[1]:

data[i, :] = temp_data[:, order]

else:

raise ValueError(f"{this} not present in provided MGXS library")

return data

def _calculate_mgxs_elem_mat(this, types, library, orders=None,

temperature=294., ce_cross_sections=None,

enrichment=None):

"""Determines the multi-group cross sections of an element or material

object.

If the data for the nuclide or macroscopic object in the library is

represented as angle-dependent data then this method will return the

geometric average cross section over all angles.

Parameters

----------

this : str or openmc.Material

Object to source data from. Elements can be input as a str

types : Iterable of str

The type of cross sections to calculate; values can either be those

in openmc.PLOT_TYPES_MGXS

library : openmc.MGXSLibrary

MGXS Library containing the data of interest

orders : Iterable of Integral, optional

The scattering order or delayed group index to use for the

corresponding entry in types. Defaults to the 0th order for scattering

and the total delayed neutron data.

temperature : float, optional

Temperature in Kelvin to plot. If not specified, a default

temperature of 294K will be plotted. Note that the nearest

temperature in the library for each nuclide will be used as opposed

to using any interpolation.

ce_cross_sections : str, optional

Location of continuous-energy cross_sections.xml file. Default is None.

This is used only for expanding the elements

enrichment : float, optional

Enrichment for U235 in weight percent. For example, input 4.95 for

4.95 weight percent enriched U. Default is None

(natural composition).

Returns

-------

data : numpy.ndarray

Cross sections calculated at the energy grid described by energy_grid

"""

if isinstance(this, openmc.Material):

if this.temperature is not None:

T = this.temperature

else:

T = temperature

# Check to see if we have nuclides/elements or a macroscopic object

if this._macroscopic is not None:

# We have macroscopics

nuclides = {this._macroscopic: this.density}

else:

# Expand elements in to nuclides with atomic densities

nuclides = this.get_nuclide_atom_densities()

# For ease of processing split out nuc and nuc_density

nuc_fraction = list(nuclides.values())

else:

T = temperature

# Expand elements in to nuclides with atomic densities

nuclides = openmc.Element(this).expand(100., 'ao', enrichment=enrichment,

cross_sections=ce_cross_sections)

# For ease of processing split out nuc and nuc_fractions

nuc_fraction = [nuclide[1] for nuclide in nuclides]

nuc_data = []

for nuclide in nuclides.items():

nuc_data.append(_calculate_mgxs_nuc_macro(nuclide[0], types, library,

orders, T))

# Combine across the nuclides

data = np.zeros((len(types), library.energy_groups.num_groups))

for line in range(len(types)):

if types[line] == 'unity':

data[line, :] = 1.

else:

for n in range(len(nuclides)):

data[line, :] += nuc_fraction[n] * nuc_data[n][line, :]

return data

openmc/openmc/plotter.py at develop · openmc-dev/openmc · GitHub

Uh oh!

Expand file tree

Breadcrumbs

plotter.py

Latest commit

History

Breadcrumbs

plotter.py

File metadata and controls