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

Open-Sn / opensn / 20185909776

12 Dec 2025 06:55PM UTC coverage: 74.333% (+0.3%) from 74.037%
20185909776

push

github

web-flow
Merge pull request #859 from wdhawkins/td_source_driver

Adding time-dependent solver and time-dependent sources

367 of 398 new or added lines in 23 files covered. (92.21%)

113 existing lines in 28 files now uncovered.

18610 of 25036 relevant lines covered (74.33%)

68947552.69 hits per line

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

0.0
/modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h
1
// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2
// SPDX-License-Identifier: MIT
3

4
#pragma once
5

6
#include "modules/problem.h"
7
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep.h"
8
#include "modules/linear_boltzmann_solvers/lbs_problem/groupset/lbs_groupset.h"
9
#include "modules/linear_boltzmann_solvers/lbs_problem/source_functions/source_flags.h"
10
#include "modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.h"
11
#include "modules/linear_boltzmann_solvers/lbs_problem/volumetric_source/volumetric_source.h"
12
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_structs.h"
13
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_view.h"
14
#include "framework/math/spatial_discretization/spatial_discretization.h"
15
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
16
#include "framework/math/linear_solver/linear_solver.h"
17
#include "framework/math/spatial_discretization/finite_element/unit_cell_matrices.h"
18
#include "framework/math/geometry.h"
19
#include <petscksp.h>
20
#include <any>
21
#include <chrono>
22

23
namespace opensn
24
{
25

26
class MPICommunicatorSet;
27
class GridFaceHistogram;
28
class AGSLinearSolver;
29
class WGSLinearSolver;
30
struct WGSContext;
31

32
/// Base class for all Linear Boltzmann Solvers.
33
class LBSProblem : public Problem
34
{
35
public:
36
  explicit LBSProblem(const std::string& name, std::shared_ptr<MeshContinuum> grid);
37

38
  /// Input parameters based construction.
39
  explicit LBSProblem(const InputParameters& params);
40

41
  LBSProblem(const LBSProblem&) = delete;
42

43
  LBSProblem& operator=(const LBSProblem&) = delete;
44

45
  ~LBSProblem() override;
46

47
  /// Returns a reference to the solver options.
48
  LBSOptions& GetOptions();
49

50
  /// Returns a constant reference to the solver options.
51
  const LBSOptions& GetOptions() const;
52

53
  void SetOptions(const InputParameters& input);
54

55
  /// Returns simulation time in seconds for time dependent problems
56
  double GetTime() const;
57

58
  /// Sets simulation time in seconds for time dependent problems
59
  void SetTime(double time);
60

61
  /// Sets dt
62
  void SetTimeStep(double dt);
63

64
  /// Returns dt
65
  double GetTimeStep() const;
66

67
  /// Sets theta for time discretization
68
  void SetTheta(double theta);
69

70
  /// Sets theta for time discretization
71
  double GetTheta() const;
72

73
  /// Is the problem time dependent
74
  bool IsTimeDependent() const;
75

76
  void SetBoundaryOptions(const InputParameters& params);
77

78
  void SetAdjoint(bool adjoint);
79

80
  GeometryType GetGeometryType() const;
81

82
  /// Returns the number of moments for the solver. This will only be non-zero after initialization.
83
  size_t GetNumMoments() const;
84

85
  unsigned int GetMaxCellDOFCount() const;
86

87
  unsigned int GetMinCellDOFCount() const;
88

89
  bool UseGPUs() const;
90

91
  /// Returns the number of groups for the solver. This will only be non-zero after initialization.
92
  size_t GetNumGroups() const;
93

94
  /// Returns the scattering order for the solver. This will only be non-zero after initialization.
95
  unsigned int GetScatteringOrder() const;
96

97
  /**
98
   * Returns the number of precursors for the solver. This will only be non-zero after
99
   * initialization.
100
   */
101
  size_t GetNumPrecursors() const;
102

103
  /**
104
   * Returns the maximum number of precursors defined on any material. This will only be non-zero
105
   * after initialization.
106
   */
107
  size_t GetMaxPrecursorsPerMaterial() const;
108

109
  const std::vector<LBSGroup>& GetGroups() const;
110

111
  std::vector<LBSGroupset>& GetGroupsets();
112

113
  const std::vector<LBSGroupset>& GetGroupsets() const;
114

115
  /// Adds a point source to the solver.
116
  void AddPointSource(std::shared_ptr<PointSource> point_source);
117

118
  /// Clears all the point sources from the solver.
119
  void ClearPointSources();
120

121
  /// Constant accessor to the list of point sources.
122
  const std::vector<std::shared_ptr<PointSource>>& GetPointSources() const;
123

124
  /// Adds a volumetric source to the solver.
125
  void AddVolumetricSource(std::shared_ptr<VolumetricSource> volumetric_source);
126

127
  /// Clears all the volumetric sources from the solver.
128
  void ClearVolumetricSources();
129

130
  /// Constant accessor to the list of volumetric sources.
131
  const std::vector<std::shared_ptr<VolumetricSource>>& GetVolumetricSources() const;
132

133
  /// Clears all the boundary conditions from the solver.
134
  void ClearBoundaries();
135

136
  size_t& GetLastRestartTime();
137

138
  /// Returns a reference to the map of material ids to XSs.
139
  const BlockID2XSMap& GetBlockID2XSMap() const;
140

141
  /// Obtains a reference to the grid.
142
  std::shared_ptr<MeshContinuum> GetGrid() const;
143

144
  /// Get pointer to carriers.
145
  void* GetCarrier(std::uint32_t idx) { return carriers_.at(idx); }
146

147
  /// Get pointer to pinners.
148
  void* GetPinner(std::uint32_t idx) { return pinners_.at(idx); }
149

150
  /// Obtains a reference to the spatial discretization.
151
  const SpatialDiscretization& GetSpatialDiscretization() const;
152

153
  /// Returns read-only access to the unit cell matrices.
154
  const std::vector<UnitCellMatrices>& GetUnitCellMatrices() const;
155

156
  /// Returns read-only access to the unit ghost cell matrices.
157
  const std::map<uint64_t, UnitCellMatrices>& GetUnitGhostCellMatrices() const;
158

159
  /// Returns a reference to the list of local cell transport views.
160
  std::vector<CellLBSView>& GetCellTransportViews();
161

162
  /// Returns a const reference to the list of local cell transport views.
163
  const std::vector<CellLBSView>& GetCellTransportViews() const;
164

165
  /// Read/Write access to the boundary preferences.
166
  std::map<uint64_t, BoundaryPreference>& GetBoundaryPreferences();
167

168
  /// Obtains a reference to the unknown manager for flux-moments.
169
  const UnknownManager& GetUnknownManager() const;
170

171
  /// Returns the local node count for the flux-moments data structures.
172
  size_t GetLocalNodeCount() const;
173

174
  /// Returns the global node count for the flux-moments data structures.
175
  size_t GetGlobalNodeCount() const;
176

177
  /// Read/write access to source moments vector.
178
  std::vector<double>& GetQMomentsLocal();
179

180
  /// Read access to source moments vector.
181
  const std::vector<double>& GetQMomentsLocal() const;
182

183
  /// Read/write access to exterior src moments vector.
184
  std::vector<double>& GetExtSrcMomentsLocal();
185

186
  /// Read access to exterior src moments vector.
187
  const std::vector<double>& GetExtSrcMomentsLocal() const;
188

189
  /// Read/write access to last updated flux vector.
190
  std::vector<double>& GetPhiOldLocal();
191

192
  /// Read access to last updated flux vector.
193
  const std::vector<double>& GetPhiOldLocal() const;
194

195
  /// Read/write access to newest updated flux vector.
196
  std::vector<double>& GetPhiNewLocal();
197

198
  /// Read access to newest updated flux vector.
199
  const std::vector<double>& GetPhiNewLocal() const;
200

201
  /// Read/write access to newest updated precursors vector.
202
  std::vector<double>& GetPrecursorsNewLocal();
203

204
  /// Read access to newest updated precursors vector.
205
  const std::vector<double>& GetPrecursorsNewLocal() const;
206

207
  /// Read/write access to the cell-wise densities.
208
  std::vector<double>& GetDensitiesLocal();
209

210
  /// Read access to the cell-wise densities.
211
  const std::vector<double>& GetDensitiesLocal() const;
212

213
  SetSourceFunction GetActiveSetSourceFunction() const;
214

215
  std::shared_ptr<AGSLinearSolver> GetAGSSolver();
216

217
  std::vector<std::shared_ptr<LinearSolver>>& GetWGSSolvers();
218

219
  WGSContext& GetWGSContext(int groupset_id);
220

221
  /**
222
   * Gets the local and global number of iterative unknowns. This normally is only the flux moments,
223
   * however, the sweep based solvers might include delayed angular fluxes in this number.
224
   */
225
  virtual std::pair<size_t, size_t> GetNumPhiIterativeUnknowns();
226

227
  /// Gets the local handle of a flux-moment based field function.
228
  size_t MapPhiFieldFunction(size_t g, size_t m) const;
229

230
  /// Returns the power generation field function, if enabled.
231
  std::shared_ptr<FieldFunctionGridBased> GetPowerFieldFunction() const;
232

233
  void Initialize() override;
234

235
  bool TriggerRestartDump() const
236
  {
237
    if (options_.write_restart_time_interval <= std::chrono::seconds(0))
238
      return false;
239

240
    auto elapsed = std::chrono::system_clock::now() - options_.last_restart_write_time;
241
    return elapsed >= options_.write_restart_time_interval;
242
  }
243

244
  void UpdateRestartWriteTime()
×
245
  {
246
    options_.last_restart_write_time = std::chrono::system_clock::now();
×
247
  }
248

249
  /// Makes a source-moments vector from scattering and fission based on the latest phi-solution.
250
  std::vector<double> MakeSourceMomentsFromPhi();
251

252
  /// Copy relevant section of phi_old to the field functions.
253
  void UpdateFieldFunctions();
254

255
  /// Sets the internal phi vector to the value in the associated field function.
256
  void SetPhiFromFieldFunctions(PhiSTLOption which_phi,
257
                                const std::vector<size_t>& m_indices,
258
                                const std::vector<size_t>& g_indices);
259

260
  /**
261
   * A method for post-processing an adjoint solution.
262
   *
263
   * @note This does nothing for diffusion-based solvers.
264
   */
265
  virtual void ReorientAdjointSolution() {};
×
266

NEW
267
  virtual void UpdatePsiOld() {};
×
268

269
protected:
270
  virtual void PrintSimHeader();
271

272
  void ComputeUnitIntegrals();
273

274
  virtual void InitializeSpatialDiscretization();
275

276
  /// Initializes parallel arrays.
277
  void InitializeParrays();
278

279
  void InitializeFieldFunctions();
280

281
  /// Initializes boundaries.
UNCOV
282
  virtual void InitializeBoundaries() {}
×
283

284
  void InitializeSolverSchemes();
285

UNCOV
286
  virtual void InitializeWGSSolvers() {};
×
UNCOV
287

×
UNCOV
288
  /// Initializes data carriers to GPUs and memory pinner.
×
UNCOV
289
  void InitializeGPUExtras();
×
UNCOV
290

×
291
  /// Reset data carriers to null and unpin memory.
292
  void ResetGPUCarriers();
293

294
  virtual void ZeroSolutions() = 0;
295

296
  LBSOptions options_;
297
  double time_ = 0.0;
298
  bool time_dependent_ = false;
299
  double theta_ = 1.0;
300
  double dt_ = 1.0;
301
  GeometryType geometry_type_ = GeometryType::INVALID;
302
  size_t num_moments_ = 0;
303
  size_t num_groups_ = 0;
304
  unsigned int scattering_order_ = 0;
305
  size_t num_precursors_ = 0;
306
  size_t max_precursors_per_material_ = 0;
307

308
  std::vector<LBSGroup> groups_;
309
  std::vector<LBSGroupset> groupsets_;
310

311
  BlockID2XSMap block_id_to_xs_map_;
312

313
  std::vector<std::shared_ptr<PointSource>> point_sources_;
314
  std::vector<std::shared_ptr<VolumetricSource>> volumetric_sources_;
315

316
  std::shared_ptr<MeshContinuum> grid_;
317
  std::shared_ptr<SpatialDiscretization> discretization_ = nullptr;
318

319
  std::vector<CellFaceNodalMapping> grid_nodal_mappings_;
320
  std::shared_ptr<MPICommunicatorSet> grid_local_comm_set_ = nullptr;
321

322
  std::vector<UnitCellMatrices> unit_cell_matrices_;
323
  std::map<uint64_t, UnitCellMatrices> unit_ghost_cell_matrices_;
324
  std::vector<CellLBSView> cell_transport_views_;
325

326
  std::map<uint64_t, BoundaryPreference> boundary_preferences_;
327

328
  UnknownManager flux_moments_uk_man_;
329

330
  unsigned int max_cell_dof_count_ = 0;
331
  unsigned int min_cell_dof_count_ = 0;
332
  uint64_t local_node_count_ = 0;
333
  uint64_t global_node_count_ = 0;
334

335
  std::vector<double> q_moments_local_, ext_src_moments_local_;
336
  std::vector<double> phi_new_local_, phi_old_local_;
337
  std::vector<double> precursor_new_local_;
338
  std::vector<double> densities_local_;
339

340
  SetSourceFunction active_set_source_function_;
341

342
  std::shared_ptr<AGSLinearSolver> ags_solver_;
343
  std::vector<std::shared_ptr<LinearSolver>> wgs_solvers_;
344

345
  std::map<std::pair<size_t, size_t>, size_t> phi_field_functions_local_map_;
346
  size_t power_gen_fieldfunc_local_handle_ = 0;
347

348
  /**
349
   * \brief Data carriers for necessary data to run the sweep on GPU.
350
   * \details These objects manage GPU memory allocation automatically, organize cross-section,
351
   * outflow, and mesh data into contiguous memory on the CPU, and handle copying it to the GPU.
352
   *
353
   * There are 3 carriers, respectively for cross sections, outflow and mesh.
354
   */
355
  std::array<void*, 3> carriers_ = {nullptr, nullptr, nullptr};
356

357
  /// Memory pinner for source moments and destination phi.
358
  std::array<void*, 2> pinners_ = {nullptr, nullptr};
359

360
  /// Flag indicating if GPU acceleration is enabled.
361
  bool use_gpus_;
362

363
private:
364
  /// Initialize groupsets
365
  void InitializeGroupsets(const InputParameters& params);
366

367
  /// Initializes materials
368
  void InitializeXSmapAndDensities(const InputParameters& params);
369
  void InitializeMaterials();
370

371
  /// Initialize sources
372
  void InitializeSources(const InputParameters& params);
373

374
  /// Initialize boundary conditions
375
  void InitializeBoundaryConditions(const InputParameters& params);
376

377
public:
378
  /// Max number of DOFs per cell that the sweep kernel on GPU can handle.
379
  static constexpr std::uint32_t max_dofs_gpu = 10;
380

381
  /// Returns the input parameters for this object.
382
  static InputParameters GetInputParameters();
383

384
  static InputParameters GetOptionsBlock();
385

386
  static InputParameters GetBoundaryOptionsBlock();
387

388
  static InputParameters GetXSMapEntryBlock();
389

390
protected:
391
  /// Checks if the current CPU is associated with any GPU.
392
  static void CheckCapableDevices();
393
};
394

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