18300593117

Committed 06 Oct 2025 10:47PM UTC coverage: 74.862% (-0.2%) from 75.031%

Build # 18300593117

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #759 from wdhawkins/performance

Sweep performance optimizations

Coverage Stats

294 of 302 new or added lines in 15 files covered. (97.35%)

334 existing lines in 80 files now uncovered.

17788 of 23761 relevant lines covered (74.86%)

61852783.95 hits per line

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

96.67

/modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/wgs_linear_solver.cc

// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
// SPDX-License-Identifier: MIT

#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/wgs_linear_solver.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/wgs_convergence_test.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_vecops.h"
#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
#include "framework/math/petsc_utils/petsc_utils.h"
#include "framework/logging/log.h"
#include "framework/utils/timer.h"
#include "framework/runtime.h"
#include <petscksp.h>
#include "caliper/cali.h"
#include <memory>
#include <iomanip>

namespace opensn
{

WGSLinearSolver::WGSLinearSolver(const std::shared_ptr<WGSContext>& gs_context_ptr)
  : PETScLinearSolver(gs_context_ptr->groupset.iterative_method, gs_context_ptr)
{
  auto& groupset = gs_context_ptr->groupset;
  auto& solver_tol_options = this->GetToleranceOptions();
  solver_tol_options.residual_absolute = groupset.residual_tolerance;
  solver_tol_options.maximum_iterations = groupset.max_iterations;
  solver_tol_options.gmres_restart_interval = groupset.gmres_restart_intvl;
}

WGSLinearSolver::~WGSLinearSolver()
{
  MatDestroy(&A_);
}

void
WGSLinearSolver::PreSetupCallback()
{
  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  gs_context_ptr->PreSetupCallback();
}

void
WGSLinearSolver::SetConvergenceTest()
{
  KSPSetConvergenceTest(ksp_, &GSConvergenceTest, nullptr, nullptr);
}

void
WGSLinearSolver::SetSystemSize()
{
  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  const auto sizes = gs_context_ptr->GetSystemSize();
  num_local_dofs_ = sizes.first;
  num_global_dofs_ = sizes.second;
}

void
WGSLinearSolver::SetSystem()
{
  if (IsSystemSet())
    return;

  x_ = CreateVector(num_local_dofs_, num_global_dofs_);

  VecSet(x_, 0.0);
  VecDuplicate(x_, &b_);

  // Create the matrix-shell
  MatCreateShell(opensn::mpi_comm,
                 num_local_dofs_,
                 num_local_dofs_,
                 num_global_dofs_,
                 num_global_dofs_,
                 &(*context_ptr_),
                 &A_);

  // Set the action-operator
  MatShellSetOperation(A_, MATOP_MULT, (void (*)())LinearSolverMatrixAction);

  // Set solver operators
  KSPSetOperators(ksp_, A_, A_);
  KSPSetUp(ksp_);
}

void
WGSLinearSolver::SetPreconditioner()
{
  if (IsSystemSet())
    return;
  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  gs_context_ptr->SetPreconditioner(ksp_);
}

void
WGSLinearSolver::PostSetupCallback()
{
  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  gs_context_ptr->PostSetupCallback();
}

void
WGSLinearSolver::PreSolveCallback()
{
  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  auto& groupset = gs_context_ptr->groupset;
  auto& do_problem = gs_context_ptr->do_problem;
  if (do_problem.GetOptions().verbose_inner_iterations)
  {
    log.Log() << "\n\n"
              << "********** Solving groupset " << groupset.id << " with "
              << this->GetIterativeMethodName() << "\n\n"
              << "Quadrature number of angles: " << groupset.quadrature->abscissae.size() << "\n"
              << "Groups " << groupset.groups.front().id << " " << groupset.groups.back().id
              << "\n\n";
  }
  gs_context_ptr->PreSolveCallback();
}

void
WGSLinearSolver::SetInitialGuess()
{
  // If the norm of the initial guess is large enough, the initial guess will be used, otherwise it
  // is assumed to be zero.

  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  auto& groupset = gs_context_ptr->groupset;
  auto& do_problem = gs_context_ptr->do_problem;

  LBSVecOps::SetGSPETScVecFromPrimarySTLvector(do_problem, groupset, x_, PhiSTLOption::PHI_OLD);

  double init_guess_norm = 0.0;
  VecNorm(x_, NORM_2, &init_guess_norm);

  if (init_guess_norm > 1.0e-10)
  {
    KSPSetInitialGuessNonzero(ksp_, PETSC_TRUE);
    if (gs_context_ptr->log_info)
      log.Log() << "Using phi_old as initial guess.";
  }
}

void
WGSLinearSolver::SetRHS()
{
  CALI_CXX_MARK_SCOPE("WGSLinearSolver::SetRHS");

  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  auto& groupset = gs_context_ptr->groupset;
  auto& do_problem = gs_context_ptr->do_problem;

  if (gs_context_ptr->log_info)
    log.Log() << program_timer.GetTimeString() << " Computing b";

  // SetSource for RHS
  saved_q_moments_local_ = do_problem.GetQMomentsLocal();

  const bool single_richardson =
    groupset.iterative_method == LinearSystemSolver::IterativeMethod::PETSC_RICHARDSON and
    tolerance_options.maximum_iterations == 1;

  if (not single_richardson)
  {
    const auto scope = gs_context_ptr->rhs_src_scope | ZERO_INCOMING_DELAYED_PSI;
    gs_context_ptr->set_source_function(
      groupset, do_problem.GetQMomentsLocal(), do_problem.GetPhiOldLocal(), scope);

    // Apply transport operator
    gs_context_ptr->ApplyInverseTransportOperator(scope);

    // Assemble PETSc vector
    LBSVecOps::SetGSPETScVecFromPrimarySTLvector(do_problem, groupset, b_, PhiSTLOption::PHI_NEW);

    // Compute RHS norm
    VecNorm(b_, NORM_2, &gs_context_ptr->rhs_norm);

    // Compute precondition RHS norm
    PC pc;
    KSPGetPC(ksp_, &pc);
    Vec temp_vec;
    VecDuplicate(b_, &temp_vec);
    PCApply(pc, b_, temp_vec);
    VecNorm(temp_vec, NORM_2, &gs_context_ptr->rhs_preconditioned_norm);
    VecDestroy(&temp_vec);
  }
  // If we have a single richardson iteration then the user probably wants
  // only a single sweep. Therefore, we are going to combine the scattering
  // source (normally included in the lhs_src_scope) into the sweep for the
  // RHS, and just suppress the kspsolve part.
  else
  {
    const auto scope = gs_context_ptr->rhs_src_scope | gs_context_ptr->lhs_src_scope;
    gs_context_ptr->set_source_function(
      groupset, do_problem.GetQMomentsLocal(), do_problem.GetPhiOldLocal(), scope);

    // Apply transport operator
    gs_context_ptr->ApplyInverseTransportOperator(scope);

    // Assemble PETSc vector
    LBSVecOps::SetGSPETScVecFromPrimarySTLvector(do_problem, groupset, x_, PhiSTLOption::PHI_NEW);

    // Compute RHS norm
    VecNorm(x_, NORM_2, &gs_context_ptr->rhs_norm);

    // Compute precondition RHS norm
    PC pc;
    KSPGetPC(ksp_, &pc);
    Vec temp_vec;
    VecDuplicate(x_, &temp_vec);
    PCApply(pc, x_, temp_vec);
    VecNorm(temp_vec, NORM_2, &gs_context_ptr->rhs_preconditioned_norm);
    VecDestroy(&temp_vec);

    SetKSPSolveSuppressionFlag(true);
  }
}

void
WGSLinearSolver::PostSolveCallback()
{
  // Get convergence reason
  if (not GetKSPSolveSuppressionFlag())
  {
    KSPConvergedReason reason;
    KSPGetConvergedReason(ksp_, &reason);
    if (reason != KSP_CONVERGED_RTOL and reason != KSP_DIVERGED_ITS)
      log.Log0Warning() << "Krylov solver failed. "
                        << "Reason: " << GetPETScConvergedReasonstring(reason);
  }

  // Copy x to local solution
  auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);
  auto& groupset = gs_context_ptr->groupset;
  auto& do_problem = gs_context_ptr->do_problem;
  LBSVecOps::SetPrimarySTLvectorFromGSPETScVec(do_problem, groupset, x_, PhiSTLOption::PHI_NEW);
  LBSVecOps::SetPrimarySTLvectorFromGSPETScVec(do_problem, groupset, x_, PhiSTLOption::PHI_OLD);

  // Restore saved q_moms
  do_problem.GetQMomentsLocal() = saved_q_moments_local_;

  // Context specific callback
  gs_context_ptr->PostSolveCallback();
}

} // namespace opensn

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/wgs_linear_solver.h"
5	#include "modules/linear_boltzmann_solvers/lbs_problem/iterative_methods/wgs_convergence_test.h"
6	#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_vecops.h"
7	#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
8	#include "framework/math/petsc_utils/petsc_utils.h"
9	#include "framework/logging/log.h"
10	#include "framework/utils/timer.h"
11	#include "framework/runtime.h"
12	#include <petscksp.h>
13	#include "caliper/cali.h"
14	#include <memory>
15	#include <iomanip>
16
17	namespace opensn
18	{
19
20	WGSLinearSolver::WGSLinearSolver(const std::shared_ptr<WGSContext>& gs_context_ptr)	315✔
21	: PETScLinearSolver(gs_context_ptr->groupset.iterative_method, gs_context_ptr)	315✔
22	{
23	auto& groupset = gs_context_ptr->groupset;	315✔
24	auto& solver_tol_options = this->GetToleranceOptions();	315✔
25	solver_tol_options.residual_absolute = groupset.residual_tolerance;	315✔
26	solver_tol_options.maximum_iterations = groupset.max_iterations;	315✔
27	solver_tol_options.gmres_restart_interval = groupset.gmres_restart_intvl;	315✔
28	}	315✔
29
30	WGSLinearSolver::~WGSLinearSolver()	311✔
31	{
32	MatDestroy(&A_);
33	}	311✔
34		311✔
35	void
36	WGSLinearSolver::PreSetupCallback()	287✔
37	{
38	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	287✔
39	gs_context_ptr->PreSetupCallback();	287✔
40	}	287✔
41
42	void
43	WGSLinearSolver::SetConvergenceTest()	287✔
44	{
45	KSPSetConvergenceTest(ksp_, &GSConvergenceTest, nullptr, nullptr);	287✔
46	}	287✔
47
48	void
49	WGSLinearSolver::SetSystemSize()	287✔
50	{
51	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	287✔
52	const auto sizes = gs_context_ptr->GetSystemSize();	287✔
53	num_local_dofs_ = sizes.first;	287✔
54	num_global_dofs_ = sizes.second;	287✔
55	}	287✔
56
57	void
58	WGSLinearSolver::SetSystem()	287✔
59	{
60	if (IsSystemSet())	287✔
UNCOV 61	return;	×
62
63	x_ = CreateVector(num_local_dofs_, num_global_dofs_);	287✔
64
65	VecSet(x_, 0.0);	287✔
66	VecDuplicate(x_, &b_);	287✔
67
68	// Create the matrix-shell
69	MatCreateShell(opensn::mpi_comm,	287✔
70	num_local_dofs_,
71	num_local_dofs_,
72	num_global_dofs_,
73	num_global_dofs_,
74	&(*context_ptr_),	287✔
75	&A_);
76
77	// Set the action-operator
78	MatShellSetOperation(A_, MATOP_MULT, (void (*)())LinearSolverMatrixAction);	287✔
79
80	// Set solver operators
81	KSPSetOperators(ksp_, A_, A_);	287✔
82	KSPSetUp(ksp_);	287✔
83	}
84
85	void
86	WGSLinearSolver::SetPreconditioner()	287✔
87	{
88	if (IsSystemSet())	287✔
89	return;	×
90	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	287✔
91	gs_context_ptr->SetPreconditioner(ksp_);	287✔
92	}	287✔
93
94	void
95	WGSLinearSolver::PostSetupCallback()	287✔
96	{
97	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	287✔
98	gs_context_ptr->PostSetupCallback();	287✔
99	}	287✔
100
101	void
102	WGSLinearSolver::PreSolveCallback()	1,148✔
103	{
104	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	1,148✔
105	auto& groupset = gs_context_ptr->groupset;	1,148✔
106	auto& do_problem = gs_context_ptr->do_problem;	1,148✔
107	if (do_problem.GetOptions().verbose_inner_iterations)	1,148✔
108	{
109	log.Log() << "\n\n"	1,948✔
110	<< "********** Solving groupset " << groupset.id << " with "	974✔
111	<< this->GetIterativeMethodName() << "\n\n"	1,948✔
112	<< "Quadrature number of angles: " << groupset.quadrature->abscissae.size() << "\n"	1,948✔
113	<< "Groups " << groupset.groups.front().id << " " << groupset.groups.back().id	974✔
114	<< "\n\n";	974✔
115	}
116	gs_context_ptr->PreSolveCallback();	1,148✔
117	}	1,148✔
118
119	void
120	WGSLinearSolver::SetInitialGuess()	1,148✔
121	{
122	// If the norm of the initial guess is large enough, the initial guess will be used, otherwise it
123	// is assumed to be zero.
124
125	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	1,148✔
126	auto& groupset = gs_context_ptr->groupset;	1,148✔
127	auto& do_problem = gs_context_ptr->do_problem;	1,148✔
128
129	LBSVecOps::SetGSPETScVecFromPrimarySTLvector(do_problem, groupset, x_, PhiSTLOption::PHI_OLD);	1,148✔
130
131	double init_guess_norm = 0.0;	1,148✔
132	VecNorm(x_, NORM_2, &init_guess_norm);	1,148✔
133
134	if (init_guess_norm > 1.0e-10)	1,148✔
135	{
136	KSPSetInitialGuessNonzero(ksp_, PETSC_TRUE);	888✔
137	if (gs_context_ptr->log_info)	888✔
138	log.Log() << "Using phi_old as initial guess.";	714✔
139	}
140	}	1,148✔
141
142	void
143	WGSLinearSolver::SetRHS()	1,148✔
144	{
145	CALI_CXX_MARK_SCOPE("WGSLinearSolver::SetRHS");	1,148✔
146
147	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	1,148✔
148	auto& groupset = gs_context_ptr->groupset;	1,148✔
149	auto& do_problem = gs_context_ptr->do_problem;	1,148✔
150
151	if (gs_context_ptr->log_info)	1,148✔
152	log.Log() << program_timer.GetTimeString() << " Computing b";	974✔
153
154	// SetSource for RHS
155	saved_q_moments_local_ = do_problem.GetQMomentsLocal();	1,148✔
156
157	const bool single_richardson =	1,148✔
158	groupset.iterative_method == LinearSystemSolver::IterativeMethod::PETSC_RICHARDSON and	1,328✔
159	tolerance_options.maximum_iterations == 1;	180✔
160
161	if (not single_richardson)	1,148✔
162	{
163	const auto scope = gs_context_ptr->rhs_src_scope \| ZERO_INCOMING_DELAYED_PSI;	1,036✔
164	gs_context_ptr->set_source_function(	2,072✔
165	groupset, do_problem.GetQMomentsLocal(), do_problem.GetPhiOldLocal(), scope);	1,036✔
166
167	// Apply transport operator
168	gs_context_ptr->ApplyInverseTransportOperator(scope);	1,036✔
169
170	// Assemble PETSc vector
171	LBSVecOps::SetGSPETScVecFromPrimarySTLvector(do_problem, groupset, b_, PhiSTLOption::PHI_NEW);	1,036✔
172
173	// Compute RHS norm
174	VecNorm(b_, NORM_2, &gs_context_ptr->rhs_norm);	1,036✔
175
176	// Compute precondition RHS norm
177	PC pc;
178	KSPGetPC(ksp_, &pc);	1,036✔
179	Vec temp_vec;
180	VecDuplicate(b_, &temp_vec);	1,036✔
181	PCApply(pc, b_, temp_vec);	1,036✔
182	VecNorm(temp_vec, NORM_2, &gs_context_ptr->rhs_preconditioned_norm);	1,036✔
183	VecDestroy(&temp_vec);	1,036✔
184	}
185	// If we have a single richardson iteration then the user probably wants
186	// only a single sweep. Therefore, we are going to combine the scattering
187	// source (normally included in the lhs_src_scope) into the sweep for the
188	// RHS, and just suppress the kspsolve part.
189	else
190	{
191	const auto scope = gs_context_ptr->rhs_src_scope \| gs_context_ptr->lhs_src_scope;	112✔
192	gs_context_ptr->set_source_function(	224✔
193	groupset, do_problem.GetQMomentsLocal(), do_problem.GetPhiOldLocal(), scope);	112✔
194
195	// Apply transport operator
196	gs_context_ptr->ApplyInverseTransportOperator(scope);	112✔
197
198	// Assemble PETSc vector
199	LBSVecOps::SetGSPETScVecFromPrimarySTLvector(do_problem, groupset, x_, PhiSTLOption::PHI_NEW);	112✔
200
201	// Compute RHS norm
202	VecNorm(x_, NORM_2, &gs_context_ptr->rhs_norm);	112✔
203
204	// Compute precondition RHS norm
205	PC pc;
206	KSPGetPC(ksp_, &pc);	112✔
207	Vec temp_vec;
208	VecDuplicate(x_, &temp_vec);	112✔
209	PCApply(pc, x_, temp_vec);	112✔
210	VecNorm(temp_vec, NORM_2, &gs_context_ptr->rhs_preconditioned_norm);	112✔
211	VecDestroy(&temp_vec);	112✔
212
213	SetKSPSolveSuppressionFlag(true);	112✔
214	}
215	}	1,148✔
216
217	void
218	WGSLinearSolver::PostSolveCallback()	1,148✔
219	{
220	// Get convergence reason
221	if (not GetKSPSolveSuppressionFlag())	1,148✔
222	{
223	KSPConvergedReason reason;
224	KSPGetConvergedReason(ksp_, &reason);	1,036✔
225	if (reason != KSP_CONVERGED_RTOL and reason != KSP_DIVERGED_ITS)	1,036✔
226	log.Log0Warning() << "Krylov solver failed. "	×
227	<< "Reason: " << GetPETScConvergedReasonstring(reason);	×
228	}
229
230	// Copy x to local solution
231	auto gs_context_ptr = std::dynamic_pointer_cast<WGSContext>(context_ptr_);	1,148✔
232	auto& groupset = gs_context_ptr->groupset;	1,148✔
233	auto& do_problem = gs_context_ptr->do_problem;	1,148✔
234	LBSVecOps::SetPrimarySTLvectorFromGSPETScVec(do_problem, groupset, x_, PhiSTLOption::PHI_NEW);	1,148✔
235	LBSVecOps::SetPrimarySTLvectorFromGSPETScVec(do_problem, groupset, x_, PhiSTLOption::PHI_OLD);	1,148✔
236
237	// Restore saved q_moms
238	do_problem.GetQMomentsLocal() = saved_q_moments_local_;	1,148✔
239
240	// Context specific callback
241	gs_context_ptr->PostSolveCallback();	1,148✔
242	}	1,148✔
243
244	} // namespace opensn

Open-Sn / opensn / 18300593117

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

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