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

Open-Sn / opensn / 21738357923

06 Feb 2026 02:45AM UTC coverage: 74.553% (+0.6%) from 73.983%
21738357923

push

github

web-flow
Merge pull request #921 from wdhawkins/update_quick_install_doc

Updating quick install instructions.

19363 of 25972 relevant lines covered (74.55%)

59700221.09 hits per line

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

83.19
/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/time_dependent_solver.h"
14
#include "modules/linear_boltzmann_solvers/solvers/transient_solver.h"
15
#include "modules/linear_boltzmann_solvers/solvers/steady_state_solver.h"
16
#include "modules/linear_boltzmann_solvers/solvers/nl_keigen_solver.h"
17
#include "modules/linear_boltzmann_solvers/solvers/pi_keigen_solver.h"
18
#include "modules/linear_boltzmann_solvers/lbs_problem/io/lbs_problem_io.h"
19
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
20
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_compute.h"
21
#include "modules/solver.h"
22
#include <pybind11/numpy.h>
23
#include <algorithm>
24
#include <cstddef>
25
#include <cstdint>
26
#include <map>
27
#include <memory>
28
#include <string>
29
#include <vector>
30

31
namespace opensn
32
{
33

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

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

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

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

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

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

117
    Wrapper of :cpp:class:`opensn::LBSProblem`.
118
    )"
119
  );
603✔
120
  lbs_problem.def(
603✔
121
    "GetScalarFieldFunctionList",
122
    [](LBSProblem& self, bool only_scalar_flux)
603✔
123
    {
124
      py::list field_function_list_per_group;
311✔
125
      for (std::size_t group = 0; group < self.GetNumGroups(); group++)
17,087✔
126
      {
127
        if (only_scalar_flux)
16,776✔
128
        {
129
          std::size_t ff_index = self.MapPhiFieldFunction(group, 0);
11,108✔
130
          field_function_list_per_group.append(self.GetFieldFunctions()[ff_index]);
11,108✔
131
        }
132
        else
133
        {
134
          py::list field_function_list_per_moment;
5,668✔
135
          for (std::size_t moment = 0; moment < self.GetNumMoments(); moment++)
22,256✔
136
          {
137
            std::size_t ff_index = self.MapPhiFieldFunction(group, moment);
16,588✔
138
            field_function_list_per_moment.append(self.GetFieldFunctions()[ff_index]);
16,588✔
139
          }
140
          field_function_list_per_group.append(field_function_list_per_moment);
5,668✔
141
        }
5,668✔
142
      }
143
      return field_function_list_per_group;
311✔
144
    },
×
145
    R"(
146
    Return field functions grouped by energy group and, optionally, by moment.
147

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

767
    Parameters
768
    ----------
769
    clear_boundary_conditions: bool, default=False
770
        If true, all current boundary conditions are deleted.
771
    boundary_conditions: List[Dict]
772
        A list of boundary condition dictionaries. Each dictionary supports:
773
          - name: str (required)
774
              Boundary name that identifies the specific boundary.
775
          - type: {'vacuum', 'isotropic', 'reflecting', 'arbitrary'} (required)
776
              Boundary type specification.
777
          - group_strength: List[float], optional
778
              Required when ``type='isotropic'``. Isotropic strength per group.
779
          - function: AngularFluxFunction, optional
780
              Required when ``type='arbitrary'``. Callable that returns incoming angular flux.
781
    )"
782
  );
783
  do_problem.def(
603✔
784
    "GetPsi",
785
    [](DiscreteOrdinatesProblem& self)
611✔
786
    {
787
      const auto& psi = self.GetPsiNewLocal();
8✔
788
      py::list psi_list;
8✔
789
      for (const auto& vec : psi)
16✔
790
      {
791
        auto array = py::array_t<double>(static_cast<py::ssize_t>(vec.size()),
8✔
792
                                         vec.data(),
793
                                         py::cast(self));
8✔
794
        psi_list.append(array);
8✔
795
      }
8✔
796
      return psi_list;
8✔
797
    },
×
798
    R"(
799
    Return psi as a list of NumPy arrays (float64), using zero-copy views into the
800
    underlying data.
801
    )"
802
  );
803
  do_problem.def(
603✔
804
    "ComputeLeakage",
805
    [](DiscreteOrdinatesProblem& self, py::list bnd_names)
634✔
806
    {
807
      auto grid = self.GetGrid();
31✔
808
      // get the supported boundaries
809
      std::map<std::string, std::uint64_t> allowed_bd_names = grid->GetBoundaryNameMap();
31✔
810
      std::map<std::uint64_t, std::string> allowed_bd_ids = grid->GetBoundaryIDMap();
31✔
811
      // get the boundaries to parse, preserving user order
812
      std::vector<std::uint64_t> bndry_ids;
31✔
813
      if (bnd_names.size() > 1)
31✔
814
      {
815
        for (py::handle name : bnd_names)
144✔
816
        {
817
          auto sname = name.cast<std::string>();
88✔
818
          bndry_ids.push_back(allowed_bd_names.at(sname));
88✔
819
        }
88✔
820
      }
821
      else
822
      {
823
        bndry_ids = self.GetGrid()->GetUniqueBoundaryIDs();
9✔
824
      }
825
      // compute the leakage
826
      std::map<std::uint64_t, std::vector<double>> leakage = ComputeLeakage(self, bndry_ids);
31✔
827
      // convert result to native Python
828
      py::dict result;
31✔
829
      for (const auto& bndry_id : bndry_ids)
125✔
830
      {
831
        const auto it = leakage.find(bndry_id);
94✔
832
        if (it == leakage.end())
94✔
833
          continue;
×
834
        // construct numpy array and copy contents
835
        const auto& grp_wise_leakage = it->second;
94✔
836
        py::array_t<double> np_vector(py::ssize_t(grp_wise_leakage.size()));
94✔
837
        auto buffer = np_vector.request();
94✔
838
        auto *np_vector_data = static_cast<double*>(buffer.ptr);
94✔
839
        std::copy(grp_wise_leakage.begin(), grp_wise_leakage.end(), np_vector_data);
94✔
840
        const std::string& name = allowed_bd_ids.at(bndry_id);
94✔
841
        result[py::str(name)] = std::move(np_vector);
188✔
842
      }
94✔
843

844
      return result;
31✔
845
    },
62✔
846
    R"(
847
    Compute leakage for the problem.
848

849
    Parameters
850
    ----------
851
    bnd_names : List[str]
852
        A list of boundary names for which leakage should be computed.
853

854
    Returns
855
    -------
856
    Dict[str, numpy.ndarray]
857
        A dictionary mapping boundary names to group-wise leakage vectors.
858
        Each array contains the outgoing angular flux (per group) integrated over
859
        the corresponding boundary surface.
860

861
    Raises
862
    ------
863
    RuntimeError
864
        If `save_angular_flux` option was not enabled during problem setup.
865

866
    ValueError
867
        If one or more boundary ids are not present on the current mesh.
868
    )",
869
    py::arg("bnd_names")
603✔
870
  );
871

872
  // discrete ordinates curvilinear problem
873
  auto do_curvilinear_problem = py::class_<DiscreteOrdinatesCurvilinearProblem,
603✔
874
                                           std::shared_ptr<DiscreteOrdinatesCurvilinearProblem>,
875
                                           DiscreteOrdinatesProblem>(
876
    slv,
877
    "DiscreteOrdinatesCurvilinearProblem",
878
    R"(
879
    Base class for discrete ordinates problems in curvilinear geometry.
880

881
    Wrapper of :cpp:class:`opensn::DiscreteOrdinatesCurvilinearProblem`.
882
    )"
883
  );
603✔
884
  do_curvilinear_problem.def(
1,206✔
885
    py::init(
603✔
886
      [](py::kwargs& params)
8✔
887
      {
888
        return DiscreteOrdinatesCurvilinearProblem::Create(kwargs_to_param_block(params));
8✔
889
      }
890
    ),
891
    R"(
892
    Construct a discrete ordinates problem for curvilinear geometry.
893

894
    Warnings
895
    --------
896
       DiscreteOrdinatesCurvilinearProblem is **experimental** and should be used with caution!
897

898
    Parameters
899
    ----------
900
    mesh : MeshContinuum
901
        The spatial mesh.
902
    coord_system : int
903
        Coordinate system to use. Must be set to 2 (cylindrical coordinates).
904
    num_groups : int
905
        The total number of energy groups.
906
    groupsets : list of dict
907
        A list of input parameter blocks, each block provides the iterative properties for a
908
        groupset. Each dictionary supports:
909
          - groups_from_to: List[int] (required)
910
              Two-entry list with the first and last group id for the groupset, e.g. ``[0, 3]``.
911
          - angular_quadrature: pyopensn.aquad.AngularQuadrature, optional
912
              Handle to an angular quadrature.
913
          - angle_aggregation_type: {'polar', 'single', 'azimuthal'}, default='polar'
914
              Angle aggregation method to use during sweeping.
915
          - angle_aggregation_num_subsets: int, default=1
916
              Number of angle subsets used for aggregation.
917
          - inner_linear_method: {'classic_richardson', 'petsc_richardson',
918
            'petsc_gmres', 'petsc_bicgstab'}, default='petsc_richardson'
919
              Iterative method used for inner linear solves.
920
          - l_abs_tol: float, default=1.0e-6
921
              Inner linear solver absolute residual tolerance.
922
          - l_max_its: int, default=200
923
              Inner linear solver maximum iterations.
924
          - gmres_restart_interval: int, default=30
925
              GMRES restart interval, if GMRES is used.
926
          - allow_cycles: bool, default=True
927
              Whether cyclic dependencies are allowed in sweeps.
928
          - apply_wgdsa: bool, default=False
929
              Enable within-group DSA for this groupset.
930
          - wgdsa_l_abs_tol: float, default=1.0e-4
931
              WGDSA linear absolute tolerance.
932
          - wgdsa_l_max_its: int, default=30
933
              WGDSA maximum iterations.
934
          - wgdsa_verbose: bool, default=False
935
              Verbose WGDSA output.
936
          - wgdsa_petsc_options: str, default=''
937
              PETSc options string for the WGDSA solver.
938
          - apply_tgdsa: bool, default=False
939
              Enable two-grid DSA for this groupset.
940
          - tgdsa_l_abs_tol: float, default=1.0e-4
941
              TGDSA linear absolute tolerance.
942
          - tgdsa_l_max_its: int, default=30
943
              TGDSA maximum iterations.
944
          - tgdsa_verbose: bool, default=False
945
              Verbose TGDSA output.
946
          - tgdsa_petsc_options: str, default=''
947
    xs_map : list of dict
948
        A list of mappings from block ids to cross-section definitions. Each dictionary supports:
949
          - block_ids: List[int] (required)
950
              Mesh block IDs to associate with the cross section.
951
          - xs: pyopensn.xs.MultiGroupXS (required)
952
              Cross-section object to assign to the specified blocks.
953
    boundary_conditions: List[Dict], default=[]
954
        A list containing tables for each boundary specification. Each dictionary supports:
955
          - name: str (required)
956
              Boundary name that identifies the specific boundary.
957
          - type: {'vacuum', 'isotropic', 'reflecting', 'arbitrary'} (required)
958
              Boundary type specification.
959
          - group_strength: List[float], optional
960
              Required when ``type='isotropic'``. Isotropic strength per group.
961
          - function: AngularFluxFunction, optional
962
              Required when ``type='arbitrary'``. Callable that returns incoming angular flux.
963
    point_sources: List[pyopensn.source.PointSource], default=[]
964
        A list of point sources.
965
    volumetric_sources: List[pyopensn.source.VolumetricSource], default=[]
966
        A list of volumetric sources.
967
    options : dict, optional
968
        A block of optional configuration parameters. Each dictionary supports the same keys as
969
        :meth:`LBSProblem.SetOptions`, including:
970
          - max_mpi_message_size: int, default=32768
971
          - restart_writes_enabled: bool, default=False
972
          - write_delayed_psi_to_restart: bool, default=True
973
          - read_restart_path: str, default=''
974
          - write_restart_path: str, default=''
975
          - write_restart_time_interval: int, default=0
976
          - use_precursors: bool, default=False
977
          - use_source_moments: bool, default=False
978
          - save_angular_flux: bool, default=False
979
          - verbose_inner_iterations: bool, default=True
980
          - verbose_outer_iterations: bool, default=True
981
          - max_ags_iterations: int, default=100
982
          - ags_tolerance: float, default=1.0e-6
983
          - ags_convergence_check: {'l2', 'pointwise'}, default='l2'
984
          - verbose_ags_iterations: bool, default=True
985
          - power_field_function_on: bool, default=False
986
          - power_default_kappa: float, default=3.20435e-11
987
          - power_normalization: float, default=-1.0
988
          - field_function_prefix_option: {'prefix', 'solver_name'}, default='prefix'
989
          - field_function_prefix: str, default=''
990
    sweep_type : str, optional
991
        The sweep type to use. Must be one of `AAH` or `CBC`. Defaults to `AAH`.
992
    )"
993
  );
994
}
603✔
995

996
// Wrap steady-state solver
997
void
998
WrapSteadyState(py::module& slv)
603✔
999
{
1000
  // clang-format off
1001
  // steady state solver
1002
  auto steady_state_solver = py::class_<SteadyStateSourceSolver, std::shared_ptr<SteadyStateSourceSolver>,
603✔
1003
                                        Solver>(
1004
    slv,
1005
    "SteadyStateSourceSolver",
1006
    R"(
1007
    Steady state solver.
1008

1009
    Wrapper of :cpp:class:`opensn::SteadyStateSourceSolver`.
1010
    )"
1011
  );
603✔
1012
  steady_state_solver.def(
1,206✔
1013
    py::init(
603✔
1014
      [](py::kwargs& params)
285✔
1015
      {
1016
        return SteadyStateSourceSolver::Create(kwargs_to_param_block(params));
285✔
1017
      }
1018
    ),
1019
    R"(
1020
    Construct a steady state solver.
1021

1022
    Parameters
1023
    ----------
1024
    pyopensn.solver.LBSProblem : LBSProblem
1025
        Existing LBSProblem instance.
1026
    )"
1027
  );
1028
  // clang-format on
1029
}
603✔
1030

1031
// Wrap transient solver
1032
void
1033
WrapTransient(py::module& slv)
603✔
1034
{
1035
  // clang-format off
1036
  auto transient_solver =
603✔
1037
    py::class_<TransientSolver, std::shared_ptr<TransientSolver>, Solver>(
1038
      slv,
1039
      "TransientSolver",
1040
      R"(
1041
      Transient solver.
1042

1043
      Wrapper of :cpp:class:`opensn::TransientSolver`.
1044
      )"
1045
    );
603✔
1046
  transient_solver.def(
1,206✔
1047
    py::init(
603✔
1048
      [](py::kwargs& params)
120✔
1049
      {
1050
        return TransientSolver::Create(kwargs_to_param_block(params));
120✔
1051
      }
1052
    ),
1053
    R"(
1054
    Construct a transient solver.
1055

1056
    Parameters
1057
    ----------
1058
    pyopensn.solver.DiscreteOrdinatesProblem : DiscreteOrdinatesProblem
1059
        Existing discrete ordinates problem instance.
1060
    dt : float, optional, default=2.0e-3
1061
        Time step size used during the simulation.
1062
    stop_time : float, optional, default=0.1
1063
        Simulation end time.
1064
    theta : float, optional, default=0.5
1065
        Time differencing scheme parameter.
1066
    initial_state : str, optional, default="existing"
1067
        Initial state for the transient solve. Allowed values: existing, zero.
1068
        In "zero" mode, the solver may initialize the problem internally if needed.
1069
    verbose : bool, optional, default=True
1070
        Enable verbose logging.
1071
    )"
1072
  );
1073
  transient_solver.def(
603✔
1074
    "SetTimeStep",
1075
    &TransientSolver::SetTimeStep,
603✔
1076
    R"(
1077
    Set the timestep size used by :meth:`Advance`.
1078

1079
    Parameters
1080
    ----------
1081
    dt : float
1082
        New timestep size.
1083
    )");
1084
  transient_solver.def(
603✔
1085
    "SetTheta",
1086
    &TransientSolver::SetTheta,
603✔
1087
    R"(
1088
    Set the theta parameter used by :meth:`Advance`.
1089

1090
    Parameters
1091
    ----------
1092
    theta : float
1093
        Theta value between 1.0e-16 and 1.
1094
    )");
1095
  transient_solver.def(
603✔
1096
    "Advance",
1097
    &TransientSolver::Advance,
603✔
1098
    R"(
1099
    Advance the solver by a single timestep.
1100

1101
    Notes
1102
    -----
1103
    You must call :meth:`Initialize` before calling :meth:`Advance` or
1104
    :meth:`Execute`.
1105
    )");
1106
  transient_solver.def(
603✔
1107
    "SetPreAdvanceCallback",
1108
    static_cast<void (TransientSolver::*)(std::function<void()>)>(
603✔
1109
      &TransientSolver::SetPreAdvanceCallback),
1110
    R"(
1111
    Register a callback that runs before each advance within :meth:`Execute`.
1112

1113
    Parameters
1114
    ----------
1115
    callback : Optional[Callable[[], None]]
1116
        Function invoked before the solver advances a timestep. Pass None to clear.
1117
        If the callback modifies the timestep, the new value is used for the
1118
        upcoming step.
1119
    )");
1120
  transient_solver.def(
603✔
1121
    "SetPreAdvanceCallback",
1122
    static_cast<void (TransientSolver::*)(std::nullptr_t)>(
603✔
1123
      &TransientSolver::SetPreAdvanceCallback),
1124
    "Clear the PreAdvance callback by passing None.");
1125
  transient_solver.def(
603✔
1126
    "SetPostAdvanceCallback",
1127
    static_cast<void (TransientSolver::*)(std::function<void()>)>(
603✔
1128
      &TransientSolver::SetPostAdvanceCallback),
1129
    R"(
1130
    Register a callback that runs after each advance within :meth:`Execute`.
1131

1132
    Parameters
1133
    ----------
1134
    callback : Optional[Callable[[], None]]
1135
        Function invoked after the solver advances a timestep. Pass None to clear.
1136
    )");
1137
  transient_solver.def(
603✔
1138
    "SetPostAdvanceCallback",
1139
    static_cast<void (TransientSolver::*)(std::nullptr_t)>(
603✔
1140
      &TransientSolver::SetPostAdvanceCallback),
1141
    "Clear the PostAdvance callback by passing None.");
1142
  slv.attr("BackwardEuler") = 1.0;
603✔
1143
  slv.attr("CrankNicolson") = 0.5;
1,206✔
1144
  // clang-format on
1145
}
603✔
1146

1147
// Wrap time-dependent solver
1148
void
1149
WrapTimeDependent(py::module& slv)
603✔
1150
{
1151
  // clang-format off
1152
  auto time_dependent_solver =
603✔
1153
    py::class_<TimeDependentSourceSolver, std::shared_ptr<TimeDependentSourceSolver>, Solver>(
1154
      slv,
1155
      "TimeDependentSourceSolver",
1156
      R"(
1157
      Time dependent solver.
1158

1159
      Wrapper of :cpp:class:`opensn::TimeDependentSourceSolver`.
1160
      )"
1161
    );
603✔
1162
  time_dependent_solver.def(
1,206✔
1163
    py::init(
603✔
1164
      [](py::kwargs& params)
44✔
1165
      {
1166
        return TimeDependentSourceSolver::Create(kwargs_to_param_block(params));
44✔
1167
      }
1168
    ),
1169
    R"(
1170
    Construct a time dependent solver.
1171

1172
    Parameters
1173
    ----------
1174
    pyopensn.solver.LBSProblem : LBSProblem
1175
        Existing LBSProblem instance.
1176
    dt : float, optional, default=1.0
1177
        Time step size used during the simulation.
1178
    stop_time : float, optional, default=1.0
1179
        Simulation end time.
1180
    )"
1181
  );
1182
  time_dependent_solver.def(
603✔
1183
    "Advance",
1184
    &TimeDependentSourceSolver::Advance,
603✔
1185
    R"(
1186
    Advance the solver by a single timestep.
1187

1188
    This method uses the configured `dt` and `theta` values and will return
1189
    immediately if the stop time has already been reached. Calling it
1190
    repeatedly allows users to write custom python time loops.
1191
    )");
1192
  time_dependent_solver.def(
603✔
1193
    "SetTimeStep",
1194
    &TimeDependentSourceSolver::SetTimeStep,
603✔
1195
    R"(
1196
    Set the timestep size used by :meth:`Advance`.
1197

1198
    Parameters
1199
    ----------
1200
    dt : float
1201
        New timestep size.
1202
    )");
1203
  time_dependent_solver.def(
603✔
1204
    "SetTheta",
1205
    &TimeDependentSourceSolver::SetTheta,
603✔
1206
    R"(
1207
    Set the theta parameter used by :meth:`Advance`.
1208

1209
    Parameters
1210
    ----------
1211
    theta : float
1212
        Theta value between 0 and 1.
1213
    )");
1214
  time_dependent_solver.def(
603✔
1215
    "SetPreAdvanceCallback",
1216
    static_cast<void (TimeDependentSourceSolver::*)(std::function<void()>)>(
603✔
1217
      &TimeDependentSourceSolver::SetPreAdvanceCallback),
1218
    R"(
1219
    Register a callback that runs before each call to :meth:`Advance`.
1220

1221
    Parameters
1222
    ----------
1223
    callback : Optional[Callable[[], None]]
1224
        Function invoked before the solver advances a timestep. Pass None to clear.
1225
        If the callback modifies the timestep, the new value is used for the
1226
        upcoming step.
1227
    )");
1228
  time_dependent_solver.def(
603✔
1229
    "SetPreAdvanceCallback",
1230
    static_cast<void (TimeDependentSourceSolver::*)(std::nullptr_t)>(
603✔
1231
      &TimeDependentSourceSolver::SetPreAdvanceCallback),
1232
    "Clear the PreAdvance callback by passing None.");
1233
  time_dependent_solver.def(
603✔
1234
    "SetPostAdvanceCallback",
1235
    static_cast<void (TimeDependentSourceSolver::*)(std::function<void()>)>(
603✔
1236
      &TimeDependentSourceSolver::SetPostAdvanceCallback),
1237
    R"(
1238
    Register a callback that runs after each call to :meth:`Advance`.
1239

1240
    Parameters
1241
    ----------
1242
    callback : Optional[Callable[[], None]]
1243
        Function invoked after the solver advances a timestep. Pass None to clear.
1244
    )");
1245
  time_dependent_solver.def(
603✔
1246
    "SetPostAdvanceCallback",
1247
    static_cast<void (TimeDependentSourceSolver::*)(std::nullptr_t)>(
603✔
1248
      &TimeDependentSourceSolver::SetPostAdvanceCallback),
1249
    "Clear the PostAdvance callback by passing None.");
1250
  slv.attr("BackwardEuler") = 1.0;
603✔
1251
  slv.attr("CrankNicolson") = 0.5;
1,206✔
1252
  // clang-format on
1253
}
603✔
1254

1255
// Wrap non-linear k-eigen solver
1256
void
1257
WrapNLKEigen(py::module& slv)
603✔
1258
{
1259
  // clang-format off
1260
  // non-linear k-eigen solver
1261
  auto non_linear_k_eigen_solver = py::class_<NonLinearKEigenSolver, std::shared_ptr<NonLinearKEigenSolver>,
603✔
1262
                                              Solver>(
1263
    slv,
1264
    "NonLinearKEigenSolver",
1265
    R"(
1266
    Non-linear k-eigenvalue solver.
1267

1268
    Wrapper of :cpp:class:`opensn::NonLinearKEigenSolver`.
1269
    )"
1270
  );
603✔
1271
  non_linear_k_eigen_solver.def(
1,206✔
1272
    py::init(
603✔
1273
      [](py::kwargs& params)
28✔
1274
      {
1275
        return NonLinearKEigenSolver::Create(kwargs_to_param_block(params));
28✔
1276
      }
1277
        ),
1278
    R"(
1279
    Construct a non-linear k-eigenvalue solver.
1280

1281
    Parameters
1282
    ----------
1283
    lbs_problem: pyopensn.solver.LBSProblem
1284
        Existing LBSProblem instance.
1285
    nl_abs_tol: float, default=1.0e-8
1286
        Non-linear absolute tolerance.
1287
    nl_rel_tol: float, default=1.0e-8
1288
        Non-linear relative tolerance.
1289
    nl_sol_tol: float, default=1.0e-50
1290
        Non-linear solution tolerance.
1291
    nl_max_its: int, default=50
1292
        Non-linear algorithm maximum iterations.
1293
    l_abs_tol: float, default=1.0e-8
1294
        Linear absolute tolerance.
1295
    l_rel_tol: float, default=1.0e-8
1296
        Linear relative tolerance.
1297
    l_div_tol: float, default=1.0e6
1298
        Linear divergence tolerance.
1299
    l_max_its: int, default=50
1300
        Linear algorithm maximum iterations.
1301
    l_gmres_restart_intvl: int, default=30
1302
        GMRES restart interval.
1303
    l_gmres_breakdown_tol: float, default=1.0e6
1304
        GMRES breakdown tolerance.
1305
    reset_phi0: bool, default=True
1306
        If true, reinitializes scalar fluxes to 1.0.
1307
    num_initial_power_iterations: int, default=0
1308
        Number of initial power iterations before the non-linear solve.
1309
    )"
1310
  );
1311
  non_linear_k_eigen_solver.def(
603✔
1312
    "GetEigenvalue",
1313
    &NonLinearKEigenSolver::GetEigenvalue,
603✔
1314
    R"(
1315
    Return the current k‑eigenvalue.
1316
    )"
1317
  );
1318
  // clang-format on
1319
}
603✔
1320

1321
// Wrap power iteration solvers
1322
void
1323
WrapPIteration(py::module& slv)
603✔
1324
{
1325
  // clang-format off
1326
  // power iteration k-eigen solver
1327
  auto pi_k_eigen_solver = py::class_<PowerIterationKEigenSolver, std::shared_ptr<PowerIterationKEigenSolver>,
603✔
1328
                                      Solver>(
1329
    slv,
1330
    "PowerIterationKEigenSolver",
1331
    R"(
1332
    Power iteration k-eigenvalue solver.
1333

1334
    Wrapper of :cpp:class:`opensn::PowerIterationKEigenSolver`.
1335
    )"
1336
  );
603✔
1337
  pi_k_eigen_solver.def(
1,206✔
1338
    py::init(
603✔
1339
      [](py::kwargs& params)
127✔
1340
      {
1341
        return PowerIterationKEigenSolver::Create(kwargs_to_param_block(params));
127✔
1342
      }
1343
    ),
1344
    R"(
1345
    Construct a power iteration k-eigen solver.
1346

1347
    Parameters
1348
    ----------
1349
    problem: pyopensn.solver.LBSProblem
1350
        Existing DiscreteOrdinatesProblem instance.
1351
    acceleration: pyopensn.solver.DiscreteOrdinatesKEigenAcceleration
1352
        Optional DiscreteOrdinatesKEigenAcceleration instance for acceleration.
1353
    max_iters: int, default = 1000
1354
        Maximum power iterations allowed.
1355
    k_tol: float, default = 1.0e-10
1356
        Tolerance on the k-eigenvalue.
1357
    reset_solution: bool, default=True
1358
        If true, initialize flux moments to 1.0.
1359
    reset_phi0: bool, default=True
1360
        If true, reinitializes scalar fluxes to 1.0.
1361
    )"
1362
  );
1363
  pi_k_eigen_solver.def(
603✔
1364
    "GetEigenvalue",
1365
    &PowerIterationKEigenSolver::GetEigenvalue,
603✔
1366
    R"(
1367
    Return the current k‑eigenvalue.
1368
    )"
1369
  );
1370
  // clang-format on
1371
}
603✔
1372

1373
// Wrap LBS solver
1374
void
1375
WrapDiscreteOrdinatesKEigenAcceleration(py::module& slv)
603✔
1376
{
1377
  // clang-format off
1378
  // discrete ordinates k-eigen acceleration base
1379
  auto acceleration = py::class_<DiscreteOrdinatesKEigenAcceleration,
603✔
1380
                                 std::shared_ptr<DiscreteOrdinatesKEigenAcceleration>>(
1381
    slv,
1382
    "DiscreteOrdinatesKEigenAcceleration",
1383
    R"(
1384
    Base class for discrete ordinates k-eigenvalue acceleration methods.
1385

1386
    Wrapper of :cpp:class:`opensn::DiscreteOrdinatesKEigenAcceleration`.
1387
    )"
1388
  );
603✔
1389
  // SCDSA acceleration
1390
  auto scdsa_acceleration = py::class_<SCDSAAcceleration,
603✔
1391
                                       std::shared_ptr<SCDSAAcceleration>,
1392
                                       DiscreteOrdinatesKEigenAcceleration>(
1393
    slv,
1394
    "SCDSAAcceleration",
1395
    R"(
1396
    Construct an SCDSA accelerator for the power iteration k-eigenvalue solver.
1397

1398
    Wrapper of :cpp:class:`opensn::SCDSAAcceleration`.
1399
    )"
1400
  );
603✔
1401
  scdsa_acceleration.def(
1,206✔
1402
    py::init(
603✔
1403
      [](py::kwargs& params)
8✔
1404
      {
1405
        return SCDSAAcceleration::Create(kwargs_to_param_block(params));
8✔
1406
      }
1407
    ),
1408
    R"(
1409
    SCDSA acceleration for the power iteration k-eigenvalue solver.
1410

1411
    Parameters
1412
    ----------
1413
    problem: pyopensn.solver.LBSProblem
1414
        Existing DiscreteOrdinatesProblem instance.
1415
    l_abs_tol: float, defauilt=1.0e-10
1416
        Absolute residual tolerance.
1417
    max_iters: int, default=100
1418
        Maximum allowable iterations.
1419
    verbose: bool, default=False
1420
        If true, enables verbose output.
1421
    petsc_options: str, default="ssss"
1422
        Additional PETSc options.
1423
    pi_max_its: int, default=50
1424
        Maximum allowable iterations for inner power iterations.
1425
    pi_k_tol: float, default=1.0e-10
1426
        k-eigenvalue tolerance for the inner power iterations.
1427
    sdm: str, default="pwld"
1428
        Spatial discretization method to use for the diffusion solver. Valid choices are:
1429
            - 'pwld' : Piecewise Linear Discontinuous
1430
            - 'pwlc' : Piecewise Linear Continuous
1431
    )"
1432
  );
1433
  // SMM acceleration
1434
  auto smm_acceleration = py::class_<SMMAcceleration,
603✔
1435
                                     std::shared_ptr<SMMAcceleration>,
1436
                                     DiscreteOrdinatesKEigenAcceleration>(
1437
    slv,
1438
    "SMMAcceleration",
1439
    R"(
1440
    Construct an SMM accelerator for the power iteration k-eigenvalue solver.
1441

1442
    Wrapper of :cpp:class:`opensn::SMMAcceleration`.
1443
    )"
1444
  );
603✔
1445
  smm_acceleration.def(
1,206✔
1446
    py::init(
603✔
1447
      [](py::kwargs& params)
4✔
1448
      {
1449
        return SMMAcceleration::Create(kwargs_to_param_block(params));
4✔
1450
      }
1451
    ),
1452
    R"(
1453
    SMM acceleration for the power iteration k-eigenvalue solver.
1454

1455
    Warnings
1456
    --------
1457
       SMM acceleration is **experimental** and should be used with caution!
1458
       SMM accleration only supports problems with isotropic scattering.
1459

1460
    Parameters
1461
    ----------
1462
    problem: pyopensn.solver.LBSProblem
1463
        Existing DiscreteOrdinatesProblem instance.
1464
    l_abs_tol: float, defauilt=1.0e-10
1465
        Absolute residual tolerance.
1466
    max_iters: int, default=100
1467
        Maximum allowable iterations.
1468
    verbose: bool, default=False
1469
        If true, enables verbose output.
1470
    petsc_options: str, default="ssss"
1471
        Additional PETSc options.
1472
    pi_max_its: int, default=50
1473
        Maximum allowable iterations for inner power iterations.
1474
    pi_k_tol: float, default=1.0e-10
1475
        k-eigenvalue tolerance for the inner power iterations.
1476
    sdm: str, default="pwld"
1477
        Spatial discretization method to use for the diffusion solver. Valid choices are:
1478
            - 'pwld' : Piecewise Linear Discontinuous
1479
            - 'pwlc' : Piecewise Linear Continuous
1480
    )"
1481
  );
1482
  // clang-format on
1483
}
603✔
1484

1485
// Wrap the solver components of OpenSn
1486
void
1487
py_solver(py::module& pyopensn)
63✔
1488
{
1489
  py::module slv = pyopensn.def_submodule("solver", "Solver module.");
63✔
1490
  WrapProblem(slv);
63✔
1491
  WrapSolver(slv);
63✔
1492
  WrapLBS(slv);
63✔
1493
  WrapSteadyState(slv);
63✔
1494
  WrapTransient(slv);
63✔
1495
  WrapTimeDependent(slv);
63✔
1496
  WrapNLKEigen(slv);
63✔
1497
  WrapDiscreteOrdinatesKEigenAcceleration(slv);
63✔
1498
  WrapPIteration(slv);
63✔
1499
}
63✔
1500

1501
} // 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