• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

Open-Sn / opensn / 21933416870

12 Feb 2026 03:12AM UTC coverage: 74.409% (-0.4%) from 74.806%
21933416870

push

github

web-flow
Merge pull request #922 from wdhawkins/named_xs

Adding ability to load user-specifed OpenMC cross sections

28 of 52 new or added lines in 4 files covered. (53.85%)

418 existing lines in 9 files now uncovered.

19839 of 26662 relevant lines covered (74.41%)

67886044.11 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

77.33
/python/lib/solver.cc
1
// SPDX-FileCopyrightText: 2025 The OpenSn Authors <https://open-sn.github.io/opensn/>
2
// SPDX-License-Identifier: MIT
3

4
#include "python/lib/py_wrappers.h"
5
#include <pybind11/functional.h>
6
#include "framework/runtime.h"
7
#include "framework/field_functions/field_function_grid_based.h"
8
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/acceleration/discrete_ordinates_keigen_acceleration.h"
9
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/acceleration/scdsa_acceleration.h"
10
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/acceleration/smm_acceleration.h"
11
#include "modules/linear_boltzmann_solvers/discrete_ordinates_curvilinear_problem/discrete_ordinates_curvilinear_problem.h"
12
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
13
#include "modules/linear_boltzmann_solvers/solvers/transient_solver.h"
14
#include "modules/linear_boltzmann_solvers/solvers/steady_state_solver.h"
15
#include "modules/linear_boltzmann_solvers/solvers/nl_keigen_solver.h"
16
#include "modules/linear_boltzmann_solvers/solvers/pi_keigen_solver.h"
17
#include "modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h"
18
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
19
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_compute.h"
20
#include "modules/solver.h"
21
#include <pybind11/numpy.h>
22
#include <algorithm>
23
#include <cstddef>
24
#include <cstdint>
25
#include <map>
26
#include <memory>
27
#include <string>
28
#include <vector>
29

30
namespace opensn
31
{
32

33
// Wrap problem
34
void
35
WrapProblem(py::module& slv)
668✔
36
{
37
  // clang-format off
38
  // problem base
39
  auto problem = py::class_<Problem, std::shared_ptr<Problem>>(
668✔
40
    slv,
41
    "Problem",
42
    R"(
43
    Base class for all problems.
44

45
    Wrapper of :cpp:class:`opensn::Problem`.
46
    )"
47
  );
668✔
48
  problem.def(
668✔
49
    "GetFieldFunctions",
50
    [](Problem& self)
668✔
51
    {
UNCOV
52
      py::list ff_list;
×
53
      std::vector<std::shared_ptr<FieldFunctionGridBased>>& cpp_ff_list = self.GetFieldFunctions();
×
54
      for (std::shared_ptr<FieldFunctionGridBased>& ff : cpp_ff_list) {
×
55
        ff_list.append(ff);
×
56
      }
UNCOV
57
      return ff_list;
×
58
    },
×
59
    R"(
60
    Get the list of field functions.
61

62
    Returns
63
    -------
64
    List[pyopensn.fieldfunc.FieldFunctionGridBased]
65
        List of grid-based field functions representing solution data such as scalar fluxes.
66
    )"
67
  );
68
  // clang-format on
69
}
668✔
70

71
// Wrap solver
72
void
73
WrapSolver(py::module& slv)
668✔
74
{
75
  // clang-format off
76
  // solver base
77
  auto solver = py::class_<Solver, std::shared_ptr<Solver> >(
668✔
78
    slv,
79
    "Solver",
80
    R"(
81
    Base class for all solvers.
82

83
    Wrapper of :cpp:class:`opensn::Solver`.
84
    )"
85
  );
668✔
86
  solver.def(
668✔
87
    "Initialize",
88
    &Solver::Initialize,
668✔
89
    "Initialize the solver."
90
  );
91
  solver.def(
668✔
92
    "Execute",
93
    &Solver::Execute,
668✔
94
    "Execute the solver."
95
  );
96
  solver.def(
668✔
97
    "Advance",
98
    &Solver::Advance,
668✔
99
    "Advance time values function."
100
  );
101
  // clang-format on
102
}
668✔
103

104
// Wrap LBS solver
105
void
106
WrapLBS(py::module& slv)
668✔
107
{
108
  // clang-format off
109
  // LBS problem
110
  auto lbs_problem = py::class_<LBSProblem, std::shared_ptr<LBSProblem>, Problem>(
668✔
111
    slv,
112
    "LBSProblem",
113
    R"(
114
    Base class for all linear Boltzmann problems.
115

116
    Wrapper of :cpp:class:`opensn::LBSProblem`.
117
    )"
118
  );
668✔
119
  lbs_problem.def(
668✔
120
    "GetScalarFieldFunctionList",
121
    [](LBSProblem& self, bool only_scalar_flux)
668✔
122
    {
123
      py::list field_function_list_per_group;
363✔
124
      for (unsigned int group = 0; group < self.GetNumGroups(); group++)
18,407✔
125
      {
126
        if (only_scalar_flux)
18,044✔
127
        {
128
          std::size_t ff_index = self.MapPhiFieldFunction(group, 0);
11,608✔
129
          field_function_list_per_group.append(self.GetFieldFunctions()[ff_index]);
11,608✔
130
        }
131
        else
132
        {
133
          py::list field_function_list_per_moment;
6,436✔
134
          for (std::size_t moment = 0; moment < self.GetNumMoments(); moment++)
23,792✔
135
          {
136
            std::size_t ff_index = self.MapPhiFieldFunction(group, moment);
17,356✔
137
            field_function_list_per_moment.append(self.GetFieldFunctions()[ff_index]);
17,356✔
138
          }
139
          field_function_list_per_group.append(field_function_list_per_moment);
6,436✔
140
        }
6,436✔
141
      }
142
      return field_function_list_per_group;
363✔
UNCOV
143
    },
×
144
    R"(
145
    Return field functions grouped by energy group and, optionally, by moment.
146

147
    Parameters
148
    ----------
149
    only_scalar_flux : bool, default=True
150
        If True, returns only the zeroth moment (scalar flux) field function for each group.
151
        The result is a flat list of field functions, one per group.
152

153
        If False, returns all moment field functions for each group.
154
        The result is a nested list where each entry corresponds to a group and contains
155
        a list of field functions for all moments (e.g., scalar flux, higher-order moments).
156

157
    Returns
158
    -------
159
    Union[List[pyopensn.fieldfunc.FieldFunctionGridBased], List[List[pyopensn.fieldfunc.FieldFunctionGridBased]]]
160
        The structure of the returned list depends on the `only_scalar_flux` flag.
161

162
    Notes
163
    -----
164
    The moment index varies more rapidly than the group index when `only_scalar_flux` is False.
165
    )",
166
    py::arg("only_scalar_flux") = true
668✔
167
  );
168
  lbs_problem.def(
668✔
169
    "GetPowerFieldFunction",
170
    &LBSProblem::GetPowerFieldFunction,
668✔
171
    R"(
172
    Returns the power generation field function, if enabled.
173
    )"
174
  );
175
  lbs_problem.def(
668✔
176
    "GetTime",
177
    &LBSProblem::GetTime,
668✔
178
    R"(
179
    Get the current simulation time in seconds.
180
    )"
181
  );
182
  lbs_problem.def(
668✔
183
    "GetTimeStep",
184
    &LBSProblem::GetTimeStep,
668✔
185
    R"(
186
    Get the current timestep size.
187
    )"
188
  );
189
  lbs_problem.def(
668✔
190
    "SetOptions",
191
    [](LBSProblem& self, py::kwargs& params)
668✔
192
    {
UNCOV
193
      InputParameters input = LBSProblem::GetOptionsBlock();
×
194
      input.AssignParameters(kwargs_to_param_block(params));
×
195
      self.SetOptions(input);
×
196
    },
×
197
    R"(
198
    Set problem options from a large list of parameters.
199

200
    Parameters
201
    ----------
202
    max_mpi_message_size: int default=32768
203
        The maximum MPI message size used during sweeps.
204
    restart_writes_enabled: bool, default=False
205
        Flag that controls writing of restart dumps.
206
    write_delayed_psi_to_restart: bool, default=True
207
        Flag that controls writing of delayed angular fluxes to restarts.
208
    read_restart_path: str, default=''
209
        Full path for reading restart dumps including file basename.
210
    write_restart_path: str, default=''
211
        Full path for writing restart dumps including file basename.
212
    write_restart_time_interval: int, default=0
213
        Time interval in seconds at which restart data is to be written.
214
    use_precursors: bool, default=False
215
        Flag for using delayed neutron precursors.
216
    use_source_moments: bool, default=False
217
        Flag for ignoring fixed sources and selectively using source moments obtained elsewhere.
218
    save_angular_flux: bool, default=False
219
        Flag indicating whether angular fluxes are to be stored or not.
220
    verbose_inner_iterations: bool, default=True
221
        Flag to control verbosity of inner iterations.
222
    verbose_outer_iterations: bool, default=True
223
        Flag to control verbosity of across-groupset iterations.
224
    max_ags_iterations: int, default=100
225
        Maximum number of across-groupset iterations.
226
    ags_tolerance: float, default=1.0e-6
227
        Across-groupset iterations tolerance.
228
    ags_convergence_check: {'l2', 'pointwise'}, default='l2'
229
        Type of convergence check for AGS iterations.
230
    verbose_ags_iterations: bool, default=True
231
        Flag to control verbosity of across-groupset iterations.
232
    power_field_function_on: bool, default=False
233
        Flag to control the creation of the power generation field function. If set to ``True``, a
234
        field function will be created with the general name ``<solver_name>_power_generation``.
235
    power_default_kappa: float, default=3.20435e-11
236
        Default ``kappa`` value (Energy released per fission) to use for power generation when cross
237
        sections do not have ``kappa`` values. Default corresponds to 200 MeV per fission.
238
    power_normalization: float, default=-1.0
239
        Power normalization factor to use. Supply a negative or zero number to turn this off.
240
    field_function_prefix_option: {'prefix', 'solver_name'}, default='prefix'
241
        Prefix option on field function names. If unset, flux field functions will be exported as
242
        ``phi_gXXX_mYYY``, where ``XXX`` is the zero-padded 3-digit group number and ``YYY`` is the
243
        zero-padded 3-digit moment.
244
    field_function_prefix: str, default=''
245
        Prefix to use on all field functions. By default, this is empty. If specified, flux moments
246
        are exported as ``prefix_phi_gXXX_mYYY``.
247
    )"
248
  );
249
  lbs_problem.def(
668✔
250
    "ComputeFissionRate",
251
    [](LBSProblem& self, const std::string& scalar_flux_iterate)
668✔
252
    {
UNCOV
253
      const std::vector<double>* phi_ptr = nullptr;
×
254
      if (scalar_flux_iterate == "old")
×
255
      {
UNCOV
256
        phi_ptr = &self.GetPhiOldLocal();
×
257
      }
UNCOV
258
      else if (scalar_flux_iterate == "new")
×
259
      {
UNCOV
260
        phi_ptr = &self.GetPhiNewLocal();
×
261
      }
262
      else
263
      {
UNCOV
264
        throw std::invalid_argument("Unknown scalar_flux_iterate value: \"" + scalar_flux_iterate + "\".");
×
265
      }
UNCOV
266
      return ComputeFissionRate(self, *phi_ptr);
×
267
    },
268
    R"(
269
    Computes the total fission rate.
270

271
    Parameters
272
    ----------
273
    scalar_flux_iterate : {'old', 'new'}
274
        Specifies which scalar flux vector to use in the calculation.
275
            - 'old': Use the previous scalar flux iterate.
276
            - 'new': Use the current scalar flux iterate.
277

278
    Returns
279
    -------
280
    float
281
        The total fission rate.
282

283
    Raises
284
    ------
285
    ValueError
286
        If `scalar_flux_iterate` is not 'old' or 'new'.
287
    )",
288
    py::arg("scalar_flux_iterate")
668✔
289
  );
290
  lbs_problem.def(
668✔
291
    "ComputeFissionProduction",
292
    [](LBSProblem& self, const std::string& scalar_flux_iterate)
1,412✔
293
    {
294
      const std::vector<double>* phi_ptr = nullptr;
744✔
295
      if (scalar_flux_iterate == "old")
744✔
296
      {
297
        phi_ptr = &self.GetPhiOldLocal();
44✔
298
      }
299
      else if (scalar_flux_iterate == "new")
700✔
300
      {
301
        phi_ptr = &self.GetPhiNewLocal();
700✔
302
      }
303
      else
304
      {
UNCOV
305
        throw std::invalid_argument("Unknown scalar_flux_iterate value: \"" + scalar_flux_iterate + "\".");
×
306
      }
307
      return ComputeFissionProduction(self, *phi_ptr);
744✔
308
    },
309
    R"(
310
    Computes the total fission production (nu*fission).
311

312
    Parameters
313
    ----------
314
    scalar_flux_iterate : {'old', 'new'}
315
        Specifies which scalar flux vector to use in the calculation.
316
            - 'old': Use the previous scalar flux iterate.
317
            - 'new': Use the current scalar flux iterate.
318

319
    Returns
320
    -------
321
    float
322
        The total fission production.
323

324
    Raises
325
    ------
326
    ValueError
327
        If `scalar_flux_iterate` is not 'old' or 'new'.
328
    )",
329
    py::arg("scalar_flux_iterate")
668✔
330
  );
331
  lbs_problem.def(
668✔
332
    "GetPhiOldLocal",
333
    [](LBSProblem& self)
716✔
334
    {
335
      return convert_vector(self.GetPhiOldLocal());
48✔
336
    },
337
    R"(
338
    Get the previous scalar flux iterate (local vector).
339

340
    Returns
341
    -------
342
    memoryview
343
        Memory view of the local old scalar flux vector.
344
    )"
345
  );
346
  lbs_problem.def(
668✔
347
    "GetPhiNewLocal",
348
    [](LBSProblem& self)
716✔
349
    {
350
      return convert_vector(self.GetPhiNewLocal());
48✔
351
    },
352
    R"(
353
    Get the current scalar flux iterate (local vector).
354

355
    Returns
356
    -------
357
    memoryview
358
        Memory view of the local new scalar flux vector.
359
    )"
360
  );
361
  lbs_problem.def(
668✔
362
    "WriteFluxMoments",
363
    [](LBSProblem& self, const std::string& file_base)
700✔
364
    {
365
      LBSSolverIO::WriteFluxMoments(self, file_base);
32✔
366
    },
367
    R"(
368
    Write flux moments to file.
369

370
    Parameters
371
    ----------
372
    file_base: str
373
        File basename.
374
    )",
375
    py::arg("file_base")
668✔
376
  );
377
  lbs_problem.def(
668✔
378
    "CreateAndWriteSourceMoments",
379
    [](LBSProblem& self, const std::string& file_base)
672✔
380
    {
381
      std::vector<double> source_moments = self.MakeSourceMomentsFromPhi();
4✔
382
      LBSSolverIO::WriteFluxMoments(self, file_base, source_moments);
4✔
383
    },
4✔
384
    R"(
385
    Write source moments from latest flux iterate to file.
386

387
    Parameters
388
    ----------
389
    file_base: str
390
        File basename.
391
    )",
392
    py::arg("file_base")
668✔
393
  );
394
  lbs_problem.def(
668✔
395
    "ReadFluxMomentsAndMakeSourceMoments",
396
    [](LBSProblem& self, const std::string& file_base, bool single_file_flag)
668✔
397
    {
UNCOV
398
      LBSSolverIO::ReadFluxMoments(self, file_base, single_file_flag, self.GetExtSrcMomentsLocal());
×
399
      log.Log() << "Making source moments from flux file.";
×
400
      std::vector<double>& temp_phi = self.GetPhiOldLocal();
×
401
      self.GetPhiOldLocal() = self.GetExtSrcMomentsLocal();
×
402
      self.GetExtSrcMomentsLocal() = self.MakeSourceMomentsFromPhi();
×
403
      self.GetPhiOldLocal() = temp_phi;
×
404
    },
×
405
    R"(
406
    Read flux moments and compute corresponding source moments.
407

408
    Parameters
409
    ----------
410
    file_base: str
411
        File basename.
412
    single_file_flag: bool
413
        True if all flux moments are in a single file.
414
    )",
415
    py::arg("file_base"),
1,336✔
416
    py::arg("single_file_flag")
668✔
417
  );
418
  lbs_problem.def(
668✔
419
    "ReadSourceMoments",
420
    [](LBSProblem& self, const std::string& file_base, bool single_file_flag)
672✔
421
    {
422
      LBSSolverIO::ReadFluxMoments(self, file_base, single_file_flag, self.GetExtSrcMomentsLocal());
4✔
423
    },
4✔
424
    R"(
425
    Read source moments from file.
426

427
    Parameters
428
    ----------
429
    file_base: str
430
        File basename.
431
    single_file_flag: bool
432
        True if all source moments are in a single file.
433
    )",
434
    py::arg("file_base"),
1,336✔
435
    py::arg("single_file_flag")
668✔
436
  );
437
  lbs_problem.def(
668✔
438
    "ReadFluxMoments",
439
    [](LBSProblem& self, const std::string& file_base, bool single_file_flag)
668✔
440
    {
UNCOV
441
      LBSSolverIO::ReadFluxMoments(self, file_base, single_file_flag);
×
442
    },
443
    R"(
444
    Read flux moment data.
445

446
    Parameters
447
    ----------
448
    file_base: str
449
        File basename.
450
    single_file_flag: bool
451
        True if all flux moments are in a single file.
452
    )",
453
    py::arg("file_base"),
1,336✔
454
    py::arg("single_file_flag")
668✔
455
  );
456
  lbs_problem.def(
668✔
457
    "WriteAngularFluxes",
458
    [](DiscreteOrdinatesProblem& self, const std::string& file_base)
672✔
459
    {
460
      LBSSolverIO::WriteAngularFluxes(self, file_base);
4✔
461
    },
462
    R"(
463
    Write angular flux data to file.
464

465
    Parameters
466
    ----------
467
    file_base: str
468
        File basename.
469
    )",
470
    py::arg("file_base")
668✔
471
  );
472
  lbs_problem.def(
668✔
473
    "ReadAngularFluxes",
474
    [](DiscreteOrdinatesProblem& self, const std::string& file_base)
672✔
475
    {
476
      LBSSolverIO::ReadAngularFluxes(self, file_base);
4✔
477
    },
478
    R"(
479
    Read angular fluxes from file.
480

481
    Parameters
482
    ----------
483
    file_base: str
484
        File basename.
485
    )",
486
    py::arg("file_base")
668✔
487
  );
488
  lbs_problem.def(
668✔
489
    "SetPointSources",
490
    [](LBSProblem& self, py::kwargs& params)
668✔
491
    {
UNCOV
492
      for (auto [key, value] : params)
×
493
      {
UNCOV
494
        auto c_key = key.cast<std::string>();
×
495
        if (c_key == "clear_point_sources")
×
496
          self.ClearPointSources();
×
497
        else if (c_key == "point_sources")
×
498
        {
UNCOV
499
          auto sources = value.cast<py::list>();
×
500
          for (auto source : sources)
×
501
          {
UNCOV
502
            self.AddPointSource(source.cast<std::shared_ptr<PointSource>>());
×
503
          }
UNCOV
504
        }
×
505
        else
UNCOV
506
          throw std::runtime_error("Invalid argument provided to SetPointSources.\n");
×
507
      }
×
508
    },
×
509
    R"(
510
    Set or clear point sources.
511

512
    Parameters
513
    ----------
514
    clear_point_sources: bool, default=False
515
        If true, all current the point sources of the problem are deleted.
516
    point_sources: List[pyopensn.source.PointSource]
517
        List of new point sources to be added to the problem.
518
    )"
519
  );
520
  lbs_problem.def(
668✔
521
    "SetVolumetricSources",
522
    [](LBSProblem& self, py::kwargs& params)
696✔
523
    {
524
      for (auto [key, value] : params)
84✔
525
      {
526
        auto c_key = key.cast<std::string>();
28✔
527
        if (c_key == "clear_volumetric_sources")
28✔
528
          self.ClearVolumetricSources();
12✔
529
        else if (c_key == "volumetric_sources")
16✔
530
        {
531
          auto sources = value.cast<py::list>();
16✔
532
          for (auto source : sources)
48✔
533
          {
534
            self.AddVolumetricSource(source.cast<std::shared_ptr<VolumetricSource>>());
32✔
535
          }
536
        }
16✔
537
        else
UNCOV
538
          throw std::runtime_error("Invalid argument provided to SetVolumetricSources.\n");
×
539
      }
28✔
540
    },
28✔
541
    R"(
542
    Set or clear volumetric sources.
543

544
    Parameters
545
    ----------
546
    clear_volumetric_sources: bool, default=False
547
        If true, all current the volumetric sources of the problem are deleted.
548
    volumetric_sources: List[pyopensn.source.VolumetricSource]
549
        List of new volumetric sources to be added to the problem.
550
    )"
551
  );
552
  lbs_problem.def(
668✔
553
    "SetXSMap",
554
    [](LBSProblem& self, py::kwargs& params)
836✔
555
    {
556
      BlockID2XSMap xs_map;
168✔
557
      for (auto [key, value] : params)
504✔
558
      {
559
        auto c_key = key.cast<std::string>();
168✔
560
        if (c_key == "xs_map")
168✔
561
        {
562
          auto xs_entries = value.cast<py::list>();
168✔
563
          for (auto entry : xs_entries)
504✔
564
          {
565
            InputParameters xs_entry_pars = LBSProblem::GetXSMapEntryBlock();
168✔
566
            xs_entry_pars.AssignParameters(pyobj_to_param_block("", entry.cast<py::dict>()));
168✔
567
            const auto& block_ids =
168✔
568
              xs_entry_pars.GetParam("block_ids").GetVectorValue<unsigned int>();
168✔
569
            auto xs = xs_entry_pars.GetSharedPtrParam<MultiGroupXS>("xs");
168✔
570
            for (const auto& block_id : block_ids)
336✔
571
              xs_map[block_id] = xs;
168✔
572
          }
168✔
573
        }
168✔
574
        else
UNCOV
575
          throw std::runtime_error("Invalid argument provided to SetXSMap.\n");
×
576
      }
168✔
577
      self.SetBlockID2XSMap(xs_map);
168✔
578
    },
168✔
579
    R"(
580
    Replace the block-id to cross-section map.
581

582
    Parameters
583
    ----------
584
    xs_map: List[Dict]
585
        A list of block-id to cross-section mapping dictionaries. Each dictionary supports:
586
          - block_ids: List[int] (required)
587
              Mesh block ids to associate with the cross section.
588
          - xs: pyopensn.xs.MultiGroupXS (required)
589
              Cross section object.
590
    )"
591
  );
592
  lbs_problem.def(
668✔
593
    "SetAdjoint",
UNCOV
594
    [](LBSProblem& self, bool adjoint)
×
595
    {
596
      self.SetAdjoint(adjoint);
12✔
597
    }
598
  );
599

600
  // discrete ordinate solver
601
  auto do_problem = py::class_<DiscreteOrdinatesProblem, std::shared_ptr<DiscreteOrdinatesProblem>,
668✔
602
                               LBSProblem>(
603
    slv,
604
    "DiscreteOrdinatesProblem",
605
    R"(
606
    Base class for discrete ordinates problems in Cartesian geometry.
607

608
    This class implements the algorithms necessary to solve a problem using the discrete ordinates method.
609
    When paired with a solver base, the result is a solver instance configured for a specific problem type
610
    (steady-state, transient, adjoint, k-eigenvalue, etc.).
611

612
    Wrapper of :cpp:class:`opensn::DiscreteOrdinatesProblem`.
613
    )"
614
  );
668✔
615
  do_problem.def(
1,336✔
616
    py::init(
668✔
617
      [](py::kwargs& params)
494✔
618
      {
619
        return DiscreteOrdinatesProblem::Create(kwargs_to_param_block(params));
494✔
620
      }
621
    ),
622
    R"(
623
    Construct a discrete ordinates problem with Cartesian geometry.
624

625
    Parameters
626
    ----------
627
    mesh : MeshContinuum
628
        The spatial mesh.
629
    num_groups : int
630
        The total number of energy groups.
631
    groupsets : List[Dict], default=[]
632
        A list of input parameter blocks, each block provides the iterative properties for a
633
        groupset. Each dictionary supports:
634
          - groups_from_to: List[int] (required)
635
              Two-entry list with the first and last group id for the groupset, e.g. ``[0, 3]``.
636
          - angular_quadrature: pyopensn.aquad.AngularQuadrature, optional
637
              Handle to an angular quadrature.
638
          - angle_aggregation_type: {'polar', 'single', 'azimuthal'}, default='polar'
639
              Angle aggregation method to use during sweeping.
640
          - angle_aggregation_num_subsets: int, default=1
641
              Number of angle subsets used for aggregation.
642
          - inner_linear_method: {'classic_richardson', 'petsc_richardson',
643
            'petsc_gmres', 'petsc_bicgstab'}, default='petsc_richardson'
644
              Iterative method used for inner linear solves.
645
          - l_abs_tol: float, default=1.0e-6
646
              Inner linear solver absolute residual tolerance.
647
          - l_max_its: int, default=200
648
              Inner linear solver maximum iterations.
649
          - gmres_restart_interval: int, default=30
650
              GMRES restart interval, if GMRES is used.
651
          - allow_cycles: bool, default=True
652
              Whether cyclic dependencies are allowed in sweeps.
653
          - apply_wgdsa: bool, default=False
654
              Enable within-group DSA for this groupset.
655
          - wgdsa_l_abs_tol: float, default=1.0e-4
656
              WGDSA linear absolute tolerance.
657
          - wgdsa_l_max_its: int, default=30
658
              WGDSA maximum iterations.
659
          - wgdsa_verbose: bool, default=False
660
              Verbose WGDSA output.
661
          - wgdsa_petsc_options: str, default=''
662
              PETSc options string for the WGDSA solver.
663
          - apply_tgdsa: bool, default=False
664
              Enable two-grid DSA for this groupset.
665
          - tgdsa_l_abs_tol: float, default=1.0e-4
666
              TGDSA linear absolute tolerance.
667
          - tgdsa_l_max_its: int, default=30
668
              TGDSA maximum iterations.
669
          - tgdsa_verbose: bool, default=False
670
              Verbose TGDSA output.
671
          - tgdsa_petsc_options: str, default=''
672
              PETSc options string for the TGDSA solver.
673
    xs_map : List[Dict], default=[]
674
        A list of mappings from block ids to cross-section definitions. Each dictionary supports:
675
          - block_ids: List[int] (required)
676
              Mesh block IDs to associate with the cross section.
677
          - xs: pyopensn.xs.MultiGroupXS (required)
678
              Cross-section object to assign to the specified blocks.
679
    boundary_conditions: List[Dict], default=[]
680
        A list containing tables for each boundary specification. Each dictionary supports:
681
          - name: str (required)
682
              Boundary name that identifies the specific boundary.
683
          - type: {'vacuum', 'isotropic', 'reflecting', 'arbitrary'} (required)
684
              Boundary type specification.
685
          - group_strength: List[float], optional
686
              Required when ``type='isotropic'``. Isotropic strength per group.
687
          - function: AngularFluxFunction, optional
688
              Required when ``type='arbitrary'``. Callable that returns incoming angular flux.
689
    point_sources: List[pyopensn.source.PointSource], default=[]
690
        A list of point sources.
691
    volumetric_sources: List[pyopensn.source.VolumetricSource], default=[]
692
        A list of volumetric sources.
693
    options : Dict, default={}
694
        A block of optional configuration parameters. Each dictionary supports the same keys as
695
        :meth:`LBSProblem.SetOptions`, including:
696
          - max_mpi_message_size: int, default=32768
697
          - restart_writes_enabled: bool, default=False
698
          - write_delayed_psi_to_restart: bool, default=True
699
          - read_restart_path: str, default=''
700
          - write_restart_path: str, default=''
701
          - write_restart_time_interval: int, default=0
702
          - use_precursors: bool, default=False
703
          - use_source_moments: bool, default=False
704
          - save_angular_flux: bool, default=False
705
          - verbose_inner_iterations: bool, default=True
706
          - verbose_outer_iterations: bool, default=True
707
          - max_ags_iterations: int, default=100
708
          - ags_tolerance: float, default=1.0e-6
709
          - ags_convergence_check: {'l2', 'pointwise'}, default='l2'
710
          - verbose_ags_iterations: bool, default=True
711
          - power_field_function_on: bool, default=False
712
          - power_default_kappa: float, default=3.20435e-11
713
          - power_normalization: float, default=-1.0
714
          - field_function_prefix_option: {'prefix', 'solver_name'}, default='prefix'
715
          - field_function_prefix: str, default=''
716
    sweep_type : str, default="AAH"
717
        The sweep type to use. Must be one of `AAH` or `CBC`. Defaults to `AAH`.
718
    use_gpus : bool, default=False
719
        A flag specifying whether GPU acceleration is used for the sweep. Currently, only ``AAH`` is
720
        supported.
721
    )"
722
  );
723
  do_problem.def(
668✔
724
    "SetOptions",
725
    [](DiscreteOrdinatesProblem& self, py::kwargs& params)
680✔
726
    {
727
      InputParameters input = DiscreteOrdinatesProblem::GetOptionsBlock();
12✔
728
      input.AssignParameters(kwargs_to_param_block(params));
12✔
729
      self.SetOptions(input);
12✔
730
    },
12✔
731
    R"(
732
    Set problem options from a large list of parameters.
733

734
    Parameters
735
    ----------
736
    adjoint: bool, default=False
737
        Flag for toggling whether the solver is in adjoint mode.
738
    )"
739
  );
740
  do_problem.def(
668✔
741
    "SetBoundaryOptions",
742
    [](DiscreteOrdinatesProblem& self, py::kwargs& params)
668✔
743
    {
UNCOV
744
      for (auto [key, value] : params)
×
745
      {
UNCOV
746
        auto c_key = key.cast<std::string>();
×
747
        if (c_key == "clear_boundary_conditions")
×
748
          self.ClearBoundaries();
×
749
        else if (c_key == "boundary_conditions")
×
750
        {
UNCOV
751
          auto boundaries = value.cast<py::list>();
×
752
          for (auto boundary : boundaries)
×
753
          {
UNCOV
754
            InputParameters input = DiscreteOrdinatesProblem::GetBoundaryOptionsBlock();
×
755
            input.AssignParameters(pyobj_to_param_block("", boundary.cast<py::dict>()));
×
756
            self.SetBoundaryOptions(input);
×
757
          }
×
758
        }
×
759
        else
UNCOV
760
          throw std::runtime_error("Invalid argument provided to SetBoundaryOptions.\n");
×
761
      }
×
762
    },
×
763
    R"(
764
    Set or clear boundary conditions.
765

766
    Parameters
767
    ----------
768
    clear_boundary_conditions: bool, default=False
769
        If true, all current boundary conditions are deleted.
770
    boundary_conditions: List[Dict]
771
        A list of boundary condition dictionaries. Each dictionary supports:
772
          - name: str (required)
773
              Boundary name that identifies the specific boundary.
774
          - type: {'vacuum', 'isotropic', 'reflecting', 'arbitrary'} (required)
775
              Boundary type specification.
776
          - group_strength: List[float], optional
777
              Required when ``type='isotropic'``. Isotropic strength per group.
778
          - function: AngularFluxFunction, optional
779
              Required when ``type='arbitrary'``. Callable that returns incoming angular flux.
780
    )"
781
  );
782
  do_problem.def(
668✔
783
    "GetPsi",
784
    [](DiscreteOrdinatesProblem& self)
676✔
785
    {
786
      const auto& psi = self.GetPsiNewLocal();
8✔
787
      py::list psi_list;
8✔
788
      for (const auto& vec : psi)
16✔
789
      {
790
        auto array = py::array_t<double>(static_cast<py::ssize_t>(vec.size()),
8✔
791
                                         vec.data(),
792
                                         py::cast(self));
8✔
793
        psi_list.append(array);
8✔
794
      }
8✔
795
      return psi_list;
8✔
UNCOV
796
    },
×
797
    R"(
798
    Return psi as a list of NumPy arrays (float64), using zero-copy views into the
799
    underlying data.
800
    )"
801
  );
802
  do_problem.def(
668✔
803
    "GetAngularFieldFunctionList",
804
    [](DiscreteOrdinatesProblem& self, py::list groups, py::list angles)
668✔
805
    {
UNCOV
806
      std::vector<unsigned int> group_ids;
×
UNCOV
807
      std::vector<size_t> angle_ids;
×
UNCOV
808
      group_ids.reserve(groups.size());
×
UNCOV
809
      angle_ids.reserve(angles.size());
×
810

UNCOV
811
      for (py::handle g : groups)
×
UNCOV
812
        group_ids.push_back(g.cast<unsigned int>());
×
UNCOV
813
      for (py::handle a : angles)
×
UNCOV
814
        angle_ids.push_back(a.cast<size_t>());
×
815

UNCOV
816
      try
×
817
      {
UNCOV
818
        auto ff_list = self.GetAngularFluxFieldFunctionList(group_ids, angle_ids);
×
UNCOV
819
        py::list out;
×
UNCOV
820
        for (const auto& ff : ff_list)
×
UNCOV
821
          out.append(ff);
×
UNCOV
822
        return out;
×
UNCOV
823
      }
×
824
      catch (const std::exception& err)
×
825
      {
826
        const std::string message = err.what();
×
UNCOV
827
        if (message.find("problem not initialized") != std::string::npos)
×
UNCOV
828
          throw std::runtime_error(
×
829
            "GetAngularFieldFunctionList requires Initialize() to be called first. "
830
            "Call solver.Initialize() (or problem.Initialize()) before requesting "
UNCOV
831
            "angular flux field functions.");
×
UNCOV
832
        throw;
×
UNCOV
833
      }
×
UNCOV
834
    },
×
835
    R"(
836
    Return field functions for selected angular flux components.
837

838
    This requires `Initialize()` to have completed; it throws if the problem is
839
    not yet initialized.
840
    Note: You must enable angular flux storage (``save_angular_flux=True``) in
841
    the problem options, or use a transient/time-dependent solver that retains
842
    angular fluxes, otherwise the field functions will remain zero.
843

844
    Example
845
    -------
846
    ```python
847
    solver.Initialize()
848
    solver.Execute()
849
    ang_ff = phys.GetAngularFieldFunctionList(groups=[0], angles=[0])
850
    ```
851

852
    For transient/time-dependent runs, call this after each timestep. Two common
853
    patterns are:
854

855
    1) Use `TransientSolver.Execute()` with a post-advance callback:
856
    ```python
857
    solver = TransientSolver(problem=phys)
858
    solver.Initialize()
859

860
    def post_advance():
861
        ang_ff = phys.GetAngularFieldFunctionList(groups=[0], angles=[0])
862
        FieldFunctionGridBased.ExportMultipleToPVTU(ang_ff, "angular_flux_t")
863

864
    solver.SetPostAdvanceCallback(post_advance)
865
    solver.Execute()
866
    ```
867

868
    2) Use a custom Python loop with `TransientSolver.Advance()`:
869
    ```python
870
    solver = TransientSolver(problem=phys)
871
    solver.Initialize()
872
    for _ in range(num_steps):
873
        solver.Advance()
874
        ang_ff = phys.GetAngularFieldFunctionList(groups=[0], angles=[0])
875
        FieldFunctionGridBased.ExportMultipleToPVTU(ang_ff, "angular_flux_t")
876
    ```
877

878
    Parameters
879
    ----------
880
    groups : List[int]
881
        Global group indices to export.
882
    angles : List[int]
883
        Angle indices within the groupset quadrature for each group.
884

885
    Returns
886
    -------
887
    List[pyopensn.fieldfunc.FieldFunctionGridBased]
888
        Field functions for the requested (group, angle) pairs.
889
    )",
890
    py::arg("groups"),
1,336✔
891
    py::arg("angles")
668✔
892
  );
893
  do_problem.def(
668✔
894
    "ComputeLeakage",
895
    [](DiscreteOrdinatesProblem& self, py::list bnd_names)
707✔
896
    {
897
      auto grid = self.GetGrid();
39✔
898
      // get the supported boundaries
899
      std::map<std::string, std::uint64_t> allowed_bd_names = grid->GetBoundaryNameMap();
39✔
900
      std::map<std::uint64_t, std::string> allowed_bd_ids = grid->GetBoundaryIDMap();
39✔
901
      const auto coord_sys = grid->GetCoordinateSystem();
39✔
902
      const auto mesh_type = grid->GetType();
39✔
903
      const auto dim = grid->GetDimension();
39✔
904
      // get the boundaries to parse, preserving user order
905
      std::vector<std::uint64_t> bndry_ids;
39✔
906
      if (bnd_names.size() > 1)
39✔
907
      {
908
        for (py::handle name : bnd_names)
184✔
909
        {
910
          auto sname = name.cast<std::string>();
112✔
911
          if (coord_sys == CoordinateSystemType::CYLINDRICAL && dim == 2)
112✔
912
          {
913
            if (sname == "xmin" || sname == "xmax" || sname == "ymin" || sname == "ymax")
24✔
UNCOV
914
              throw std::runtime_error("ComputeLeakage: Boundary name '" + sname +
×
915
                                       "' is invalid for cylindrical orthogonal meshes. "
UNCOV
916
                                       "Use rmin, rmax, zmin, zmax.");
×
917

918
            if (mesh_type == MeshType::ORTHOGONAL)
24✔
919
            {
920
              if (sname == "rmin") sname = "xmin";
24✔
921
              else if (sname == "rmax") sname = "xmax";
24✔
922
              else if (sname == "zmin") sname = "ymin";
16✔
923
              else if (sname == "zmax") sname = "ymax";
8✔
924
            }
925
          }
926
          bndry_ids.push_back(allowed_bd_names.at(sname));
112✔
927
        }
112✔
928
      }
929
      else
930
      {
931
        bndry_ids = self.GetGrid()->GetUniqueBoundaryIDs();
9✔
932
      }
933
      // compute the leakage
934
      std::map<std::uint64_t, std::vector<double>> leakage = ComputeLeakage(self, bndry_ids);
39✔
935
      const bool rz_ortho = (coord_sys == CoordinateSystemType::CYLINDRICAL &&
78✔
936
                             mesh_type == MeshType::ORTHOGONAL && dim == 2);
39✔
937

938
      auto to_rz_name = [](const std::string& name)
63✔
939
      {
940
        if (name == "xmin") return std::string("rmin");
24✔
941
        if (name == "xmax") return std::string("rmax");
24✔
942
        if (name == "ymin") return std::string("zmin");
16✔
943
        if (name == "ymax") return std::string("zmax");
8✔
UNCOV
944
        return name;
×
945
      };
946

947
      // convert result to native Python
948
      py::dict result;
39✔
949
      for (const auto& bndry_id : bndry_ids)
157✔
950
      {
951
        const auto it = leakage.find(bndry_id);
118✔
952
        if (it == leakage.end())
118✔
UNCOV
953
          continue;
×
954
        // construct numpy array and copy contents
955
        const auto& grp_wise_leakage = it->second;
118✔
956
        py::array_t<double> np_vector(py::ssize_t(grp_wise_leakage.size()));
118✔
957
        auto buffer = np_vector.request();
118✔
958
        auto *np_vector_data = static_cast<double*>(buffer.ptr);
118✔
959
        std::copy(grp_wise_leakage.begin(), grp_wise_leakage.end(), np_vector_data);
118✔
960
        std::string name = allowed_bd_ids.at(bndry_id);
118✔
961
        if (rz_ortho)
118✔
962
          name = to_rz_name(name);
24✔
963
        result[py::str(name)] = std::move(np_vector);
236✔
964
      }
118✔
965

966
      return result;
39✔
967
    },
78✔
968
    R"(
969
    Compute leakage for the problem.
970

971
    Parameters
972
    ----------
973
    bnd_names : List[str]
974
        A list of boundary names for which leakage should be computed.
975

976
    Returns
977
    -------
978
    Dict[str, numpy.ndarray]
979
        A dictionary mapping boundary names to group-wise leakage vectors.
980
        Each array contains the outgoing angular flux (per group) integrated over
981
        the corresponding boundary surface.
982

983
    Raises
984
    ------
985
    RuntimeError
986
        If `save_angular_flux` option was not enabled during problem setup.
987

988
    ValueError
989
        If one or more boundary ids are not present on the current mesh.
990
    )",
991
    py::arg("bnd_names")
668✔
992
  );
993

994
  // discrete ordinates curvilinear problem
995
  auto do_curvilinear_problem = py::class_<DiscreteOrdinatesCurvilinearProblem,
668✔
996
                                           std::shared_ptr<DiscreteOrdinatesCurvilinearProblem>,
997
                                           DiscreteOrdinatesProblem>(
998
    slv,
999
    "DiscreteOrdinatesCurvilinearProblem",
1000
    R"(
1001
    Base class for discrete ordinates problems in curvilinear geometry.
1002

1003
    Wrapper of :cpp:class:`opensn::DiscreteOrdinatesCurvilinearProblem`.
1004
    )"
1005
  );
668✔
1006
  do_curvilinear_problem.def(
1,336✔
1007
    py::init(
668✔
1008
      [](py::kwargs& params)
80✔
1009
      {
1010
        return DiscreteOrdinatesCurvilinearProblem::Create(kwargs_to_param_block(params));
80✔
1011
      }
1012
    ),
1013
    R"(
1014
    Construct a discrete ordinates problem for curvilinear geometry.
1015

1016
    Warnings
1017
    --------
1018
       DiscreteOrdinatesCurvilinearProblem is **experimental** and should be used with caution!
1019

1020
    Parameters
1021
    ----------
1022
    mesh : MeshContinuum
1023
        The spatial mesh.
1024
    coord_system : int
1025
        Coordinate system to use. Must be set to 2 (cylindrical coordinates).
1026
    num_groups : int
1027
        The total number of energy groups.
1028
    groupsets : list of dict
1029
        A list of input parameter blocks, each block provides the iterative properties for a
1030
        groupset. Each dictionary supports:
1031
          - groups_from_to: List[int] (required)
1032
              Two-entry list with the first and last group id for the groupset, e.g. ``[0, 3]``.
1033
          - angular_quadrature: pyopensn.aquad.AngularQuadrature, optional
1034
              Handle to an angular quadrature.
1035
          - angle_aggregation_type: {'polar', 'single', 'azimuthal'}, default='polar'
1036
              Angle aggregation method to use during sweeping.
1037
          - angle_aggregation_num_subsets: int, default=1
1038
              Number of angle subsets used for aggregation.
1039
          - inner_linear_method: {'classic_richardson', 'petsc_richardson',
1040
            'petsc_gmres', 'petsc_bicgstab'}, default='petsc_richardson'
1041
              Iterative method used for inner linear solves.
1042
          - l_abs_tol: float, default=1.0e-6
1043
              Inner linear solver absolute residual tolerance.
1044
          - l_max_its: int, default=200
1045
              Inner linear solver maximum iterations.
1046
          - gmres_restart_interval: int, default=30
1047
              GMRES restart interval, if GMRES is used.
1048
          - allow_cycles: bool, default=True
1049
              Whether cyclic dependencies are allowed in sweeps.
1050
          - apply_wgdsa: bool, default=False
1051
              Enable within-group DSA for this groupset.
1052
          - wgdsa_l_abs_tol: float, default=1.0e-4
1053
              WGDSA linear absolute tolerance.
1054
          - wgdsa_l_max_its: int, default=30
1055
              WGDSA maximum iterations.
1056
          - wgdsa_verbose: bool, default=False
1057
              Verbose WGDSA output.
1058
          - wgdsa_petsc_options: str, default=''
1059
              PETSc options string for the WGDSA solver.
1060
          - apply_tgdsa: bool, default=False
1061
              Enable two-grid DSA for this groupset.
1062
          - tgdsa_l_abs_tol: float, default=1.0e-4
1063
              TGDSA linear absolute tolerance.
1064
          - tgdsa_l_max_its: int, default=30
1065
              TGDSA maximum iterations.
1066
          - tgdsa_verbose: bool, default=False
1067
              Verbose TGDSA output.
1068
          - tgdsa_petsc_options: str, default=''
1069
    xs_map : list of dict
1070
        A list of mappings from block ids to cross-section definitions. Each dictionary supports:
1071
          - block_ids: List[int] (required)
1072
              Mesh block IDs to associate with the cross section.
1073
          - xs: pyopensn.xs.MultiGroupXS (required)
1074
              Cross-section object to assign to the specified blocks.
1075
    boundary_conditions: List[Dict], default=[]
1076
        A list containing tables for each boundary specification. Each dictionary supports:
1077
          - name: str (required)
1078
              Boundary name that identifies the specific boundary.
1079
          - type: {'vacuum', 'isotropic', 'reflecting', 'arbitrary'} (required)
1080
              Boundary type specification.
1081
          - group_strength: List[float], optional
1082
              Required when ``type='isotropic'``. Isotropic strength per group.
1083
          - function: AngularFluxFunction, optional
1084
              Required when ``type='arbitrary'``. Callable that returns incoming angular flux.
1085
    point_sources: List[pyopensn.source.PointSource], default=[]
1086
        A list of point sources.
1087
    volumetric_sources: List[pyopensn.source.VolumetricSource], default=[]
1088
        A list of volumetric sources.
1089
    options : dict, optional
1090
        A block of optional configuration parameters. Each dictionary supports the same keys as
1091
        :meth:`LBSProblem.SetOptions`, including:
1092
          - max_mpi_message_size: int, default=32768
1093
          - restart_writes_enabled: bool, default=False
1094
          - write_delayed_psi_to_restart: bool, default=True
1095
          - read_restart_path: str, default=''
1096
          - write_restart_path: str, default=''
1097
          - write_restart_time_interval: int, default=0
1098
          - use_precursors: bool, default=False
1099
          - use_source_moments: bool, default=False
1100
          - save_angular_flux: bool, default=False
1101
          - verbose_inner_iterations: bool, default=True
1102
          - verbose_outer_iterations: bool, default=True
1103
          - max_ags_iterations: int, default=100
1104
          - ags_tolerance: float, default=1.0e-6
1105
          - ags_convergence_check: {'l2', 'pointwise'}, default='l2'
1106
          - verbose_ags_iterations: bool, default=True
1107
          - power_field_function_on: bool, default=False
1108
          - power_default_kappa: float, default=3.20435e-11
1109
          - power_normalization: float, default=-1.0
1110
          - field_function_prefix_option: {'prefix', 'solver_name'}, default='prefix'
1111
          - field_function_prefix: str, default=''
1112
    sweep_type : str, optional
1113
        The sweep type to use. Must be one of `AAH` or `CBC`. Defaults to `AAH`.
1114
    )"
1115
  );
1116
}
668✔
1117

1118
// Wrap steady-state solver
1119
void
1120
WrapSteadyState(py::module& slv)
668✔
1121
{
1122
  // clang-format off
1123
  // steady state solver
1124
  auto steady_state_solver = py::class_<SteadyStateSourceSolver, std::shared_ptr<SteadyStateSourceSolver>,
668✔
1125
                                        Solver>(
1126
    slv,
1127
    "SteadyStateSourceSolver",
1128
    R"(
1129
    Steady state solver.
1130

1131
    Wrapper of :cpp:class:`opensn::SteadyStateSourceSolver`.
1132
    )"
1133
  );
668✔
1134
  steady_state_solver.def(
1,336✔
1135
    py::init(
668✔
1136
      [](py::kwargs& params)
353✔
1137
      {
1138
        return SteadyStateSourceSolver::Create(kwargs_to_param_block(params));
353✔
1139
      }
1140
    ),
1141
    R"(
1142
    Construct a steady state solver.
1143

1144
    Parameters
1145
    ----------
1146
    pyopensn.solver.LBSProblem : LBSProblem
1147
        Existing LBSProblem instance.
1148
    )"
1149
  );
1150
  // clang-format on
1151
}
668✔
1152

1153
// Wrap transient solver
1154
void
1155
WrapTransient(py::module& slv)
668✔
1156
{
1157
  // clang-format off
1158
  auto transient_solver =
668✔
1159
    py::class_<TransientSolver, std::shared_ptr<TransientSolver>, Solver>(
1160
      slv,
1161
      "TransientSolver",
1162
      R"(
1163
      Transient solver.
1164

1165
      Wrapper of :cpp:class:`opensn::TransientSolver`.
1166
      )"
1167
    );
668✔
1168
  transient_solver.def(
1,336✔
1169
    py::init(
668✔
1170
      [](py::kwargs& params)
156✔
1171
      {
1172
        return TransientSolver::Create(kwargs_to_param_block(params));
156✔
1173
      }
1174
    ),
1175
    R"(
1176
    Construct a transient solver.
1177

1178
    Parameters
1179
    ----------
1180
    pyopensn.solver.DiscreteOrdinatesProblem : DiscreteOrdinatesProblem
1181
        Existing discrete ordinates problem instance.
1182
    dt : float, optional, default=2.0e-3
1183
        Time step size used during the simulation.
1184
    stop_time : float, optional, default=0.1
1185
        Simulation end time.
1186
    theta : float, optional, default=0.5
1187
        Time differencing scheme parameter.
1188
    initial_state : str, optional, default="existing"
1189
        Initial state for the transient solve. Allowed values: existing, zero.
1190
        In "zero" mode, the solver may initialize the problem internally if needed.
1191
    verbose : bool, optional, default=True
1192
        Enable verbose logging.
1193
    )"
1194
  );
1195
  transient_solver.def(
668✔
1196
    "SetTimeStep",
1197
    &TransientSolver::SetTimeStep,
668✔
1198
    R"(
1199
    Set the timestep size used by :meth:`Advance`.
1200

1201
    Parameters
1202
    ----------
1203
    dt : float
1204
        New timestep size.
1205
    )");
1206
  transient_solver.def(
668✔
1207
    "SetTheta",
1208
    &TransientSolver::SetTheta,
668✔
1209
    R"(
1210
    Set the theta parameter used by :meth:`Advance`.
1211

1212
    Parameters
1213
    ----------
1214
    theta : float
1215
        Theta value between 1.0e-16 and 1.
1216
    )");
1217
  transient_solver.def(
668✔
1218
    "Advance",
1219
    &TransientSolver::Advance,
668✔
1220
    R"(
1221
    Advance the solver by a single timestep.
1222

1223
    Notes
1224
    -----
1225
    You must call :meth:`Initialize` before calling :meth:`Advance` or
1226
    :meth:`Execute`.
1227
    )");
1228
  transient_solver.def(
668✔
1229
    "SetPreAdvanceCallback",
1230
    static_cast<void (TransientSolver::*)(std::function<void()>)>(
668✔
1231
      &TransientSolver::SetPreAdvanceCallback),
1232
    R"(
1233
    Register a callback that runs before each advance within :meth:`Execute`.
1234

1235
    Parameters
1236
    ----------
1237
    callback : Optional[Callable[[], None]]
1238
        Function invoked before the solver advances a timestep. Pass None to clear.
1239
        If the callback modifies the timestep, the new value is used for the
1240
        upcoming step.
1241
    )");
1242
  transient_solver.def(
668✔
1243
    "SetPreAdvanceCallback",
1244
    static_cast<void (TransientSolver::*)(std::nullptr_t)>(
668✔
1245
      &TransientSolver::SetPreAdvanceCallback),
1246
    "Clear the PreAdvance callback by passing None.");
1247
  transient_solver.def(
668✔
1248
    "SetPostAdvanceCallback",
1249
    static_cast<void (TransientSolver::*)(std::function<void()>)>(
668✔
1250
      &TransientSolver::SetPostAdvanceCallback),
1251
    R"(
1252
    Register a callback that runs after each advance within :meth:`Execute`.
1253

1254
    Parameters
1255
    ----------
1256
    callback : Optional[Callable[[], None]]
1257
        Function invoked after the solver advances a timestep. Pass None to clear.
1258
    )");
1259
  transient_solver.def(
668✔
1260
    "SetPostAdvanceCallback",
1261
    static_cast<void (TransientSolver::*)(std::nullptr_t)>(
668✔
1262
      &TransientSolver::SetPostAdvanceCallback),
1263
    "Clear the PostAdvance callback by passing None.");
1264
  slv.attr("BackwardEuler") = 1.0;
668✔
1265
  slv.attr("CrankNicolson") = 0.5;
1,336✔
1266
  // clang-format on
1267
}
668✔
1268

1269
// Wrap non-linear k-eigen solver
1270
void
1271
WrapNLKEigen(py::module& slv)
668✔
1272
{
1273
  // clang-format off
1274
  // non-linear k-eigen solver
1275
  auto non_linear_k_eigen_solver = py::class_<NonLinearKEigenSolver, std::shared_ptr<NonLinearKEigenSolver>,
668✔
1276
                                              Solver>(
1277
    slv,
1278
    "NonLinearKEigenSolver",
1279
    R"(
1280
    Non-linear k-eigenvalue solver.
1281

1282
    Wrapper of :cpp:class:`opensn::NonLinearKEigenSolver`.
1283
    )"
1284
  );
668✔
1285
  non_linear_k_eigen_solver.def(
1,336✔
1286
    py::init(
668✔
1287
      [](py::kwargs& params)
28✔
1288
      {
1289
        return NonLinearKEigenSolver::Create(kwargs_to_param_block(params));
28✔
1290
      }
1291
        ),
1292
    R"(
1293
    Construct a non-linear k-eigenvalue solver.
1294

1295
    Parameters
1296
    ----------
1297
    lbs_problem: pyopensn.solver.LBSProblem
1298
        Existing LBSProblem instance.
1299
    nl_abs_tol: float, default=1.0e-8
1300
        Non-linear absolute tolerance.
1301
    nl_rel_tol: float, default=1.0e-8
1302
        Non-linear relative tolerance.
1303
    nl_sol_tol: float, default=1.0e-50
1304
        Non-linear solution tolerance.
1305
    nl_max_its: int, default=50
1306
        Non-linear algorithm maximum iterations.
1307
    l_abs_tol: float, default=1.0e-8
1308
        Linear absolute tolerance.
1309
    l_rel_tol: float, default=1.0e-8
1310
        Linear relative tolerance.
1311
    l_div_tol: float, default=1.0e6
1312
        Linear divergence tolerance.
1313
    l_max_its: int, default=50
1314
        Linear algorithm maximum iterations.
1315
    l_gmres_restart_intvl: int, default=30
1316
        GMRES restart interval.
1317
    l_gmres_breakdown_tol: float, default=1.0e6
1318
        GMRES breakdown tolerance.
1319
    reset_phi0: bool, default=True
1320
        If true, reinitializes scalar fluxes to 1.0.
1321
    num_initial_power_iterations: int, default=0
1322
        Number of initial power iterations before the non-linear solve.
1323
    )"
1324
  );
1325
  non_linear_k_eigen_solver.def(
668✔
1326
    "GetEigenvalue",
1327
    &NonLinearKEigenSolver::GetEigenvalue,
668✔
1328
    R"(
1329
    Return the current k‑eigenvalue.
1330
    )"
1331
  );
1332
  // clang-format on
1333
}
668✔
1334

1335
// Wrap power iteration solvers
1336
void
1337
WrapPIteration(py::module& slv)
668✔
1338
{
1339
  // clang-format off
1340
  // power iteration k-eigen solver
1341
  auto pi_k_eigen_solver = py::class_<PowerIterationKEigenSolver, std::shared_ptr<PowerIterationKEigenSolver>,
668✔
1342
                                      Solver>(
1343
    slv,
1344
    "PowerIterationKEigenSolver",
1345
    R"(
1346
    Power iteration k-eigenvalue solver.
1347

1348
    Wrapper of :cpp:class:`opensn::PowerIterationKEigenSolver`.
1349
    )"
1350
  );
668✔
1351
  pi_k_eigen_solver.def(
1,336✔
1352
    py::init(
668✔
1353
      [](py::kwargs& params)
147✔
1354
      {
1355
        return PowerIterationKEigenSolver::Create(kwargs_to_param_block(params));
147✔
1356
      }
1357
    ),
1358
    R"(
1359
    Construct a power iteration k-eigen solver.
1360

1361
    Parameters
1362
    ----------
1363
    problem: pyopensn.solver.LBSProblem
1364
        Existing DiscreteOrdinatesProblem instance.
1365
    acceleration: pyopensn.solver.DiscreteOrdinatesKEigenAcceleration
1366
        Optional DiscreteOrdinatesKEigenAcceleration instance for acceleration.
1367
    max_iters: int, default = 1000
1368
        Maximum power iterations allowed.
1369
    k_tol: float, default = 1.0e-10
1370
        Tolerance on the k-eigenvalue.
1371
    reset_solution: bool, default=True
1372
        If true, initialize flux moments to 1.0.
1373
    reset_phi0: bool, default=True
1374
        If true, reinitializes scalar fluxes to 1.0.
1375
    )"
1376
  );
1377
  pi_k_eigen_solver.def(
668✔
1378
    "GetEigenvalue",
1379
    &PowerIterationKEigenSolver::GetEigenvalue,
668✔
1380
    R"(
1381
    Return the current k‑eigenvalue.
1382
    )"
1383
  );
1384
  // clang-format on
1385
}
668✔
1386

1387
// Wrap LBS solver
1388
void
1389
WrapDiscreteOrdinatesKEigenAcceleration(py::module& slv)
668✔
1390
{
1391
  // clang-format off
1392
  // discrete ordinates k-eigen acceleration base
1393
  auto acceleration = py::class_<DiscreteOrdinatesKEigenAcceleration,
668✔
1394
                                 std::shared_ptr<DiscreteOrdinatesKEigenAcceleration>>(
1395
    slv,
1396
    "DiscreteOrdinatesKEigenAcceleration",
1397
    R"(
1398
    Base class for discrete ordinates k-eigenvalue acceleration methods.
1399

1400
    Wrapper of :cpp:class:`opensn::DiscreteOrdinatesKEigenAcceleration`.
1401
    )"
1402
  );
668✔
1403
  // SCDSA acceleration
1404
  auto scdsa_acceleration = py::class_<SCDSAAcceleration,
668✔
1405
                                       std::shared_ptr<SCDSAAcceleration>,
1406
                                       DiscreteOrdinatesKEigenAcceleration>(
1407
    slv,
1408
    "SCDSAAcceleration",
1409
    R"(
1410
    Construct an SCDSA accelerator for the power iteration k-eigenvalue solver.
1411

1412
    Wrapper of :cpp:class:`opensn::SCDSAAcceleration`.
1413
    )"
1414
  );
668✔
1415
  scdsa_acceleration.def(
1,336✔
1416
    py::init(
668✔
1417
      [](py::kwargs& params)
8✔
1418
      {
1419
        return SCDSAAcceleration::Create(kwargs_to_param_block(params));
8✔
1420
      }
1421
    ),
1422
    R"(
1423
    SCDSA acceleration for the power iteration k-eigenvalue solver.
1424

1425
    Parameters
1426
    ----------
1427
    problem: pyopensn.solver.LBSProblem
1428
        Existing DiscreteOrdinatesProblem instance.
1429
    l_abs_tol: float, defauilt=1.0e-10
1430
        Absolute residual tolerance.
1431
    max_iters: int, default=100
1432
        Maximum allowable iterations.
1433
    verbose: bool, default=False
1434
        If true, enables verbose output.
1435
    petsc_options: str, default="ssss"
1436
        Additional PETSc options.
1437
    pi_max_its: int, default=50
1438
        Maximum allowable iterations for inner power iterations.
1439
    pi_k_tol: float, default=1.0e-10
1440
        k-eigenvalue tolerance for the inner power iterations.
1441
    sdm: str, default="pwld"
1442
        Spatial discretization method to use for the diffusion solver. Valid choices are:
1443
            - 'pwld' : Piecewise Linear Discontinuous
1444
            - 'pwlc' : Piecewise Linear Continuous
1445
    )"
1446
  );
1447
  // SMM acceleration
1448
  auto smm_acceleration = py::class_<SMMAcceleration,
668✔
1449
                                     std::shared_ptr<SMMAcceleration>,
1450
                                     DiscreteOrdinatesKEigenAcceleration>(
1451
    slv,
1452
    "SMMAcceleration",
1453
    R"(
1454
    Construct an SMM accelerator for the power iteration k-eigenvalue solver.
1455

1456
    Wrapper of :cpp:class:`opensn::SMMAcceleration`.
1457
    )"
1458
  );
668✔
1459
  smm_acceleration.def(
1,336✔
1460
    py::init(
668✔
1461
      [](py::kwargs& params)
4✔
1462
      {
1463
        return SMMAcceleration::Create(kwargs_to_param_block(params));
4✔
1464
      }
1465
    ),
1466
    R"(
1467
    SMM acceleration for the power iteration k-eigenvalue solver.
1468

1469
    Warnings
1470
    --------
1471
       SMM acceleration is **experimental** and should be used with caution!
1472
       SMM accleration only supports problems with isotropic scattering.
1473

1474
    Parameters
1475
    ----------
1476
    problem: pyopensn.solver.LBSProblem
1477
        Existing DiscreteOrdinatesProblem instance.
1478
    l_abs_tol: float, defauilt=1.0e-10
1479
        Absolute residual tolerance.
1480
    max_iters: int, default=100
1481
        Maximum allowable iterations.
1482
    verbose: bool, default=False
1483
        If true, enables verbose output.
1484
    petsc_options: str, default="ssss"
1485
        Additional PETSc options.
1486
    pi_max_its: int, default=50
1487
        Maximum allowable iterations for inner power iterations.
1488
    pi_k_tol: float, default=1.0e-10
1489
        k-eigenvalue tolerance for the inner power iterations.
1490
    sdm: str, default="pwld"
1491
        Spatial discretization method to use for the diffusion solver. Valid choices are:
1492
            - 'pwld' : Piecewise Linear Discontinuous
1493
            - 'pwlc' : Piecewise Linear Continuous
1494
    )"
1495
  );
1496
  // clang-format on
1497
}
668✔
1498

1499
// Wrap the solver components of OpenSn
1500
void
1501
py_solver(py::module& pyopensn)
63✔
1502
{
1503
  py::module slv = pyopensn.def_submodule("solver", "Solver module.");
63✔
1504
  WrapProblem(slv);
63✔
1505
  WrapSolver(slv);
63✔
1506
  WrapLBS(slv);
63✔
1507
  WrapSteadyState(slv);
63✔
1508
  WrapTransient(slv);
63✔
1509
  WrapNLKEigen(slv);
63✔
1510
  WrapDiscreteOrdinatesKEigenAcceleration(slv);
63✔
1511
  WrapPIteration(slv);
63✔
1512
}
63✔
1513

1514
} // namespace opensn
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc