27259186575

Committed 10 Jun 2026 06:11AM UTC coverage: 77.473% (+0.2%) from 77.306%

Build # 27259186575

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1084 from wdhawkins/cmfd_tutorial

Adding CMFD tutorial.

Coverage Stats

24792 of 32001 relevant lines covered (77.47%)

82088528.83 hits per line

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

82.11

/modules/diffusion/diffusion.cc

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

#include "modules/diffusion/diffusion.h"
#include "framework/logging/log_format.h"
#include "framework/math/spatial_discretization/spatial_discretization.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/math/petsc_utils/petsc_utils.h"
#include "framework/logging/log.h"
#include "framework/runtime.h"
#include "framework/utils/timer.h"
#include "framework/utils/error.h"
#include <sstream>

namespace opensn
{

namespace
{

void
LogDiffusionSolveStart(const std::string& name, Vec rhs)
{
  double rhs_norm = 0.0;
  OpenSnPETScCall(VecNorm(rhs, NORM_2, &rhs_norm));

  std::ostringstream out;
  out << name << " solve";
  AppendNumericField(out, "rhs_norm", rhs_norm, Scientific(6));
  log.Log() << no_wrap << program_timer.GetTimeString() << " " << out.str();
}

void
LogDiffusionSolveFinal(const std::string& name, KSP ksp, Vec x)
{
  double solution_norm = 0.0;
  OpenSnPETScCall(VecNorm(x, NORM_2, &solution_norm));

  KSPConvergedReason reason = KSP_CONVERGED_ITERATING;
  OpenSnPETScCall(KSPGetConvergedReason(ksp, &reason));
  PetscInt iterations = 0;
  OpenSnPETScCall(KSPGetIterationNumber(ksp, &iterations));

  const auto status = KSPReasonToPETScSolverStatus(reason);

  std::ostringstream out;
  out << name << " final, status = " << PETScSolverStatusName(status);
  AppendNumericField(out, "iterations", static_cast<std::size_t>(iterations));
  AppendNumericField(out, "solution_norm", solution_norm, Scientific(6));
  if (log.GetVerbosity() >= 2)
    out << ", detail = " << GetPETScConvergedReasonstring(reason);
  log.Log() << no_wrap << program_timer.GetTimeString() << " " << out.str();
}

} // namespace

DiffusionSolver::DiffusionSolver(std::string name,
                                 const opensn::SpatialDiscretization& sdm,
                                 const UnknownManager& uk_man,
                                 std::map<uint64_t, BoundaryCondition> bcs,
                                 MatID2XSMap map_mat_id_2_xs,
                                 const std::vector<UnitCellMatrices>& unit_cell_matrices,
                                 const bool suppress_bcs,
                                 const bool requires_ghosts,
                                 const bool verbose)
  : name_(std::move(name)),
    grid_(sdm.GetGrid()),
    sdm_(sdm),
    uk_man_(uk_man),
    bcs_(std::move(bcs)),
    mat_id_2_xs_map_(std::move(map_mat_id_2_xs)),
    unit_cell_matrices_(unit_cell_matrices),
    num_local_dofs_(static_cast<PetscInt>(sdm_.GetNumLocalDOFs(uk_man_))),
    num_global_dofs_(static_cast<PetscInt>(sdm_.GetNumGlobalDOFs(uk_man_))),
    A_(nullptr),
    rhs_(nullptr),
    ksp_(nullptr),
    requires_ghosts_(requires_ghosts),
    suppress_bcs_(suppress_bcs)
{
  options.verbose = verbose;
}

DiffusionSolver::~DiffusionSolver()
{
  OpenSnPETScCall(MatDestroy(&A_));
  OpenSnPETScCall(VecDestroy(&rhs_));
  OpenSnPETScCall(VecDestroy(&x_));
  OpenSnPETScCall(KSPDestroy(&ksp_));
}

std::string
DiffusionSolver::GetName() const
{
  return name_;
}

const Vec&
DiffusionSolver::GetRHS() const
{
  return rhs_;
}

const UnknownManager&
DiffusionSolver::GetUnknownStructure() const
{
  return uk_man_;
}

const SpatialDiscretization&
DiffusionSolver::GetSpatialDiscretization() const
{
  return sdm_;
}

std::pair<std::uint64_t, std::uint64_t>
DiffusionSolver::GetNumPhiIterativeUnknowns()
{
  return {sdm_.GetNumLocalDOFs(uk_man_), sdm_.GetNumGlobalDOFs(uk_man_)};
}

void
DiffusionSolver::AddToRHS(const std::vector<double>& values)
{
  const auto num_local_dofs = sdm_.GetNumLocalDOFs(uk_man_);
  if (num_local_dofs != values.size())
    throw std::invalid_argument("Vector size mismatch.");

  PetscScalar* rhs_ptr = nullptr;
  OpenSnPETScCall(VecGetArray(rhs_, &rhs_ptr));
  for (std::uint64_t i = 0; i < num_local_dofs; ++i)
    rhs_ptr[i] += values[i];
  OpenSnPETScCall(VecRestoreArray(rhs_, &rhs_ptr));
}

void
DiffusionSolver::AddToMatrix(const std::vector<PetscInt>& rows,
                             const std::vector<PetscInt>& cols,
                             const std::vector<double>& vals)
{
  if (rows.size() != cols.size() or rows.size() != vals.size())
    throw std::invalid_argument("The number of row entries, column entries, and value "
                                "entries do not agree.");
  for (int i = 0; i < vals.size(); ++i)
    OpenSnPETScCall(MatSetValue(A_, rows[i], cols[i], vals[i], ADD_VALUES));
  OpenSnPETScCall(MatAssemblyBegin(A_, MAT_FLUSH_ASSEMBLY));
  OpenSnPETScCall(MatAssemblyEnd(A_, MAT_FLUSH_ASSEMBLY));
}

void
DiffusionSolver::Initialize()
{
  if (options.verbose)
    log.Log() << name_ << ": Initializing PETSc items";

  if (options.verbose)
    log.Log() << name_ << ": Global number of DOFs=" << num_global_dofs_;

  opensn::mpi_comm.barrier();
  log.Log() << "Sparsity pattern";
  opensn::mpi_comm.barrier();
  // Create Matrix
  std::vector<int64_t> nodal_nnz_in_diag;
  std::vector<int64_t> nodal_nnz_off_diag;
  sdm_.BuildSparsityPattern(nodal_nnz_in_diag, nodal_nnz_off_diag, uk_man_);
  opensn::mpi_comm.barrier();
  log.Log() << "Done Sparsity pattern";
  opensn::mpi_comm.barrier();
  A_ = CreateSquareMatrix(num_local_dofs_, num_global_dofs_);
  InitMatrixSparsity(A_, nodal_nnz_in_diag, nodal_nnz_off_diag);
  opensn::mpi_comm.barrier();
  log.Log() << "Done matrix creation";
  opensn::mpi_comm.barrier();

  // Create RHS
  if (not requires_ghosts_)
    rhs_ = CreateVector(num_local_dofs_, num_global_dofs_);
  else
  {
    auto ghost_ids = sdm_.GetGhostDOFIndices(uk_man_);
    std::vector<int64_t> ghids(ghost_ids.begin(), ghost_ids.end());
    rhs_ = CreateVectorWithGhosts(num_local_dofs_,
                                  num_global_dofs_,
                                  static_cast<PetscInt>(sdm_.GetNumGhostDOFs(uk_man_)),
                                  ghids);
  }

  opensn::mpi_comm.barrier();
  log.Log() << "Done vector creation";
  opensn::mpi_comm.barrier();

  OpenSnPETScCall(VecDuplicate(rhs_, &x_));

  // Create KSP
  OpenSnPETScCall(KSPCreate(opensn::mpi_comm, &ksp_));
  OpenSnPETScCall(KSPSetOptionsPrefix(ksp_, name_.c_str()));

  OpenSnPETScCall(
    KSPSetTolerances(ksp_, 1.0e-50, options.residual_tolerance, 1.0e50, options.max_iters));

  // Set Pre-conditioner
  PC pc = nullptr;
  OpenSnPETScCall(KSPGetPC(ksp_, &pc));
  const bool petsc_options_policy = options.solver_policy == "petsc_options";
  const bool direct_policy = options.solver_policy == "direct" or
                             (options.solver_policy == "auto" and opensn::mpi_comm.size() == 1 and
                              num_global_dofs_ <= options.direct_solve_threshold);
  OpenSnInvalidArgumentIf(options.solver_policy != "auto" and options.solver_policy != "direct" and
                            options.solver_policy != "iterative" and not petsc_options_policy,
                          "Unsupported diffusion solver policy \"" + options.solver_policy +
                            "\". Valid values are auto, direct, iterative, and petsc_options.");

  if (not petsc_options_policy and direct_policy)
  {
    OpenSnPETScCall(KSPSetType(ksp_, KSPPREONLY));
    OpenSnPETScCall(PCSetType(pc, PCLU));
  }
  else if (not petsc_options_policy)
  {
    OpenSnPETScCall(KSPSetType(ksp_, KSPCG));
    OpenSnPETScCall(PCSetType(pc, PCHYPRE));

    OpenSnPETScCall(PCHYPRESetType(pc, "boomeramg"));
    std::vector<std::string> pc_options = {"pc_hypre_boomeramg_agg_nl 1",
                                           "pc_hypre_boomeramg_P_max 4",
                                           "pc_hypre_boomeramg_grid_sweeps_coarse 1",
                                           "pc_hypre_boomeramg_max_levels 25",
                                           "pc_hypre_boomeramg_relax_type_all symmetric-SOR/Jacobi",
                                           "pc_hypre_boomeramg_coarsen_type HMIS",
                                           "pc_hypre_boomeramg_interp_type ext+i"};

    if (grid_->GetDimension() == 2)
      pc_options.emplace_back("pc_hypre_boomeramg_strong_threshold 0.6");
    else if (grid_->GetDimension() == 3)
      pc_options.emplace_back("pc_hypre_boomeramg_strong_threshold 0.8");

    for (const auto& option : pc_options)
      OpenSnPETScCall(PetscOptionsInsertString(nullptr, ("-" + name_ + option).c_str()));
  }

  OpenSnPETScCall(PetscOptionsInsertString(nullptr, options.additional_options_string.c_str()));

  OpenSnPETScCall(PCSetFromOptions(pc));
  OpenSnPETScCall(KSPSetFromOptions(ksp_));
}

void
DiffusionSolver::Solve(std::vector<double>& solution, bool use_initial_guess)
{
  const std::string fname = "acceleration::DiffusionMIPSolver::Solve";
  OpenSnPETScCall(VecSet(x_, 0.0));

  if (not use_initial_guess)
    OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_FALSE));
  else
    OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_TRUE));

  OpenSnPETScCall(KSPSetTolerances(
    ksp_, options.residual_tolerance, options.residual_tolerance, 1.0e50, options.max_iters));

  if (options.perform_symmetry_check)
  {
    PetscBool symmetry = PETSC_FALSE;
    OpenSnPETScCall(MatIsSymmetric(A_, 1.0e-6, &symmetry));
    if (symmetry == PETSC_FALSE)
      throw std::logic_error(fname + ":Symmetry check failed");
  }

  if (options.verbose)
  {
    OpenSnPETScCall(KSPMonitorSet(ksp_, &KSPMonitorRelativeToRHS, nullptr, nullptr));
    LogDiffusionSolveStart(name_, rhs_);
  }

  if (use_initial_guess)
  {
    PetscScalar* x_raw = nullptr;
    OpenSnPETScCall(VecGetArray(x_, &x_raw));
    size_t k = 0;
    for (const auto& value : solution)
      x_raw[k++] = value;
    OpenSnPETScCall(VecRestoreArray(x_, &x_raw));
  }

  // Solve
  OpenSnPETScCall(KSPSolve(ksp_, rhs_, x_));

  // Print convergence info
  if (options.verbose)
    LogDiffusionSolveFinal(name_, ksp_, x_);

  // Transfer petsc solution to vector
  if (requires_ghosts_)
  {
    CommunicateGhostEntries(x_);
    sdm_.LocalizePETScVectorWithGhosts(x_, solution, uk_man_);
  }
  else
    sdm_.LocalizePETScVector(x_, solution, uk_man_);
}

void
DiffusionSolver::Solve(Vec petsc_solution, bool use_initial_guess)
{
  const std::string fname = "acceleration::DiffusionMIPSolver::Solve";
  OpenSnPETScCall(VecSet(x_, 0.0));

  if (not use_initial_guess)
    OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_FALSE));
  else
    OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_TRUE));

  OpenSnPETScCall(KSPSetTolerances(
    ksp_, options.residual_tolerance, options.residual_tolerance, 1.0e50, options.max_iters));

  if (options.perform_symmetry_check)
  {
    PetscBool symmetry = PETSC_FALSE;
    OpenSnPETScCall(MatIsSymmetric(A_, 1.0e-6, &symmetry));
    if (symmetry == PETSC_FALSE)
      throw std::logic_error(fname + ":Symmetry check failed");
  }

  if (options.verbose)
  {
    OpenSnPETScCall(KSPMonitorSet(ksp_, &KSPMonitorRelativeToRHS, nullptr, nullptr));
    LogDiffusionSolveStart(name_, rhs_);
  }

  if (use_initial_guess)
  {
    OpenSnPETScCall(VecCopy(petsc_solution, x_));
  }

  // Solve
  OpenSnPETScCall(KSPSolve(ksp_, rhs_, x_));

  // Print convergence info
  if (options.verbose)
    LogDiffusionSolveFinal(name_, ksp_, x_);

  // Transfer petsc solution to vector
  OpenSnPETScCall(VecCopy(x_, petsc_solution));
}

} // namespace opensn

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "modules/diffusion/diffusion.h"
5	#include "framework/logging/log_format.h"
6	#include "framework/math/spatial_discretization/spatial_discretization.h"
7	#include "framework/mesh/mesh_continuum/mesh_continuum.h"
8	#include "framework/math/petsc_utils/petsc_utils.h"
9	#include "framework/logging/log.h"
10	#include "framework/runtime.h"
11	#include "framework/utils/timer.h"
12	#include "framework/utils/error.h"
13	#include <sstream>
14
15	namespace opensn
16	{
17
18	namespace
19	{
20
21	void
22	LogDiffusionSolveStart(const std::string& name, Vec rhs)	142✔
23	{
24	double rhs_norm = 0.0;	142✔
25	OpenSnPETScCall(VecNorm(rhs, NORM_2, &rhs_norm));	142✔
26
27	std::ostringstream out;	142✔
28	out << name << " solve";	142✔
29	AppendNumericField(out, "rhs_norm", rhs_norm, Scientific(6));	142✔
30	log.Log() << no_wrap << program_timer.GetTimeString() << " " << out.str();	568✔
31	}	142✔
32
33	void
34	LogDiffusionSolveFinal(const std::string& name, KSP ksp, Vec x)	142✔
35	{
36	double solution_norm = 0.0;	142✔
37	OpenSnPETScCall(VecNorm(x, NORM_2, &solution_norm));	142✔
38
39	KSPConvergedReason reason = KSP_CONVERGED_ITERATING;	142✔
40	OpenSnPETScCall(KSPGetConvergedReason(ksp, &reason));	142✔
41	PetscInt iterations = 0;	142✔
42	OpenSnPETScCall(KSPGetIterationNumber(ksp, &iterations));	142✔
43
44	const auto status = KSPReasonToPETScSolverStatus(reason);	142✔
45
46	std::ostringstream out;	142✔
47	out << name << " final, status = " << PETScSolverStatusName(status);	142✔
48	AppendNumericField(out, "iterations", static_cast<std::size_t>(iterations));	142✔
49	AppendNumericField(out, "solution_norm", solution_norm, Scientific(6));	142✔
50	if (log.GetVerbosity() >= 2)	142✔
51	out << ", detail = " << GetPETScConvergedReasonstring(reason);	×
52	log.Log() << no_wrap << program_timer.GetTimeString() << " " << out.str();	568✔
53	}	142✔
54
55	} // namespace
56
57	DiffusionSolver::DiffusionSolver(std::string name,	874✔
58	const opensn::SpatialDiscretization& sdm,
59	const UnknownManager& uk_man,
60	std::map<uint64_t, BoundaryCondition> bcs,
61	MatID2XSMap map_mat_id_2_xs,
62	const std::vector<UnitCellMatrices>& unit_cell_matrices,
63	const bool suppress_bcs,
64	const bool requires_ghosts,
65	const bool verbose)	874✔
66	: name_(std::move(name)),	874✔
67	grid_(sdm.GetGrid()),	874✔
68	sdm_(sdm),	874✔
69	uk_man_(uk_man),	874✔
70	bcs_(std::move(bcs)),	874✔
71	mat_id_2_xs_map_(std::move(map_mat_id_2_xs)),	874✔
72	unit_cell_matrices_(unit_cell_matrices),	874✔
73	num_local_dofs_(static_cast<PetscInt>(sdm_.GetNumLocalDOFs(uk_man_))),	874✔
74	num_global_dofs_(static_cast<PetscInt>(sdm_.GetNumGlobalDOFs(uk_man_))),	874✔
75	A_(nullptr),	874✔
76	rhs_(nullptr),	874✔
77	ksp_(nullptr),	874✔
78	requires_ghosts_(requires_ghosts),	874✔
79	suppress_bcs_(suppress_bcs)	874✔
80	{
81	options.verbose = verbose;	874✔
82	}	874✔
83
84	DiffusionSolver::~DiffusionSolver()	5,052✔
85	{
86	OpenSnPETScCall(MatDestroy(&A_));
87	OpenSnPETScCall(VecDestroy(&rhs_));
88	OpenSnPETScCall(VecDestroy(&x_));
89	OpenSnPETScCall(KSPDestroy(&ksp_));
90	}	1,684✔
91		842✔
92	std::string	842✔
93	DiffusionSolver::GetName() const	842✔
94	{	842✔
95	return name_;	×
96	}	×
97
98	const Vec&
99	DiffusionSolver::GetRHS() const	×
100	{
101	return rhs_;	×
102	}
103
104	const UnknownManager&
105	DiffusionSolver::GetUnknownStructure() const	38,672✔
106	{
107	return uk_man_;	38,672✔
108	}
109
110	const SpatialDiscretization&
111	DiffusionSolver::GetSpatialDiscretization() const	13,672✔
112	{
113	return sdm_;	13,672✔
114	}
115
116	std::pair<std::uint64_t, std::uint64_t>
117	DiffusionSolver::GetNumPhiIterativeUnknowns()	×
118	{
119	return {sdm_.GetNumLocalDOFs(uk_man_), sdm_.GetNumGlobalDOFs(uk_man_)};	×
120	}
121
122	void
123	DiffusionSolver::AddToRHS(const std::vector<double>& values)	128✔
124	{
125	const auto num_local_dofs = sdm_.GetNumLocalDOFs(uk_man_);	128✔
126	if (num_local_dofs != values.size())	128✔
127	throw std::invalid_argument("Vector size mismatch.");	×
128
129	PetscScalar* rhs_ptr = nullptr;	128✔
130	OpenSnPETScCall(VecGetArray(rhs_, &rhs_ptr));	128✔
131	for (std::uint64_t i = 0; i < num_local_dofs; ++i)	2,688,128✔
132	rhs_ptr[i] += values[i];	2,688,000✔
133	OpenSnPETScCall(VecRestoreArray(rhs_, &rhs_ptr));	128✔
134	}	128✔
135
136	void
137	DiffusionSolver::AddToMatrix(const std::vector<PetscInt>& rows,	4✔
138	const std::vector<PetscInt>& cols,
139	const std::vector<double>& vals)
140	{
141	if (rows.size() != cols.size() or rows.size() != vals.size())	4✔
142	throw std::invalid_argument("The number of row entries, column entries, and value "	×
143	"entries do not agree.");	×
144	for (int i = 0; i < vals.size(); ++i)	4✔
145	OpenSnPETScCall(MatSetValue(A_, rows[i], cols[i], vals[i], ADD_VALUES));	×
146	OpenSnPETScCall(MatAssemblyBegin(A_, MAT_FLUSH_ASSEMBLY));	4✔
147	OpenSnPETScCall(MatAssemblyEnd(A_, MAT_FLUSH_ASSEMBLY));	4✔
148	}	4✔
149
150	void
151	DiffusionSolver::Initialize()	874✔
152	{
153	if (options.verbose)	874✔
154	log.Log() << name_ << ": Initializing PETSc items";	22✔
155
156	if (options.verbose)	874✔
157	log.Log() << name_ << ": Global number of DOFs=" << num_global_dofs_;	22✔
158
159	opensn::mpi_comm.barrier();	874✔
160	log.Log() << "Sparsity pattern";	1,748✔
161	opensn::mpi_comm.barrier();	874✔
162	// Create Matrix
163	std::vector<int64_t> nodal_nnz_in_diag;	874✔
164	std::vector<int64_t> nodal_nnz_off_diag;	874✔
165	sdm_.BuildSparsityPattern(nodal_nnz_in_diag, nodal_nnz_off_diag, uk_man_);	874✔
166	opensn::mpi_comm.barrier();	874✔
167	log.Log() << "Done Sparsity pattern";	1,748✔
168	opensn::mpi_comm.barrier();	874✔
169	A_ = CreateSquareMatrix(num_local_dofs_, num_global_dofs_);	874✔
170	InitMatrixSparsity(A_, nodal_nnz_in_diag, nodal_nnz_off_diag);	874✔
171	opensn::mpi_comm.barrier();	874✔
172	log.Log() << "Done matrix creation";	1,748✔
173	opensn::mpi_comm.barrier();	874✔
174
175	// Create RHS
176	if (not requires_ghosts_)	874✔
177	rhs_ = CreateVector(num_local_dofs_, num_global_dofs_);	869✔
178	else
179	{
180	auto ghost_ids = sdm_.GetGhostDOFIndices(uk_man_);	5✔
181	std::vector<int64_t> ghids(ghost_ids.begin(), ghost_ids.end());	5✔
182	rhs_ = CreateVectorWithGhosts(num_local_dofs_,	5✔
183	num_global_dofs_,	5✔
184	static_cast<PetscInt>(sdm_.GetNumGhostDOFs(uk_man_)),	5✔
185	ghids);
186	}	5✔
187
188	opensn::mpi_comm.barrier();	874✔
189	log.Log() << "Done vector creation";	1,748✔
190	opensn::mpi_comm.barrier();	874✔
191
192	OpenSnPETScCall(VecDuplicate(rhs_, &x_));	874✔
193
194	// Create KSP
195	OpenSnPETScCall(KSPCreate(opensn::mpi_comm, &ksp_));	874✔
196	OpenSnPETScCall(KSPSetOptionsPrefix(ksp_, name_.c_str()));	874✔
197
198	OpenSnPETScCall(	874✔
199	KSPSetTolerances(ksp_, 1.0e-50, options.residual_tolerance, 1.0e50, options.max_iters));
200
201	// Set Pre-conditioner
202	PC pc = nullptr;	874✔
203	OpenSnPETScCall(KSPGetPC(ksp_, &pc));	874✔
204	const bool petsc_options_policy = options.solver_policy == "petsc_options";	874✔
205	const bool direct_policy = options.solver_policy == "direct" or	874✔
206	(options.solver_policy == "auto" and opensn::mpi_comm.size() == 1 and	1,748✔
207	num_global_dofs_ <= options.direct_solve_threshold);	10✔
208	OpenSnInvalidArgumentIf(options.solver_policy != "auto" and options.solver_policy != "direct" and	874✔
209	options.solver_policy != "iterative" and not petsc_options_policy,
210	"Unsupported diffusion solver policy \"" + options.solver_policy +
211	"\". Valid values are auto, direct, iterative, and petsc_options.");
212
213	if (not petsc_options_policy and direct_policy)	874✔
214	{
215	OpenSnPETScCall(KSPSetType(ksp_, KSPPREONLY));	9✔
216	OpenSnPETScCall(PCSetType(pc, PCLU));	9✔
217	}
218	else if (not petsc_options_policy)	865✔
219	{
220	OpenSnPETScCall(KSPSetType(ksp_, KSPCG));	865✔
221	OpenSnPETScCall(PCSetType(pc, PCHYPRE));	865✔
222
223	OpenSnPETScCall(PCHYPRESetType(pc, "boomeramg"));	865✔
224	std::vector<std::string> pc_options = {"pc_hypre_boomeramg_agg_nl 1",	865✔
225	"pc_hypre_boomeramg_P_max 4",
226	"pc_hypre_boomeramg_grid_sweeps_coarse 1",
227	"pc_hypre_boomeramg_max_levels 25",
228	"pc_hypre_boomeramg_relax_type_all symmetric-SOR/Jacobi",
229	"pc_hypre_boomeramg_coarsen_type HMIS",
230	"pc_hypre_boomeramg_interp_type ext+i"};	865✔
231
232	if (grid_->GetDimension() == 2)	865✔
233	pc_options.emplace_back("pc_hypre_boomeramg_strong_threshold 0.6");	837✔
234	else if (grid_->GetDimension() == 3)	28✔
235	pc_options.emplace_back("pc_hypre_boomeramg_strong_threshold 0.8");	16✔
236
237	for (const auto& option : pc_options)	7,773✔
238	OpenSnPETScCall(PetscOptionsInsertString(nullptr, ("-" + name_ + option).c_str()));	6,908✔
239	}	865✔
240
241	OpenSnPETScCall(PetscOptionsInsertString(nullptr, options.additional_options_string.c_str()));	874✔
242
243	OpenSnPETScCall(PCSetFromOptions(pc));	874✔
244	OpenSnPETScCall(KSPSetFromOptions(ksp_));	874✔
245	}	874✔
246
247	void
248	DiffusionSolver::Solve(std::vector<double>& solution, bool use_initial_guess)	23,174✔
249	{
250	const std::string fname = "acceleration::DiffusionMIPSolver::Solve";	23,174✔
251	OpenSnPETScCall(VecSet(x_, 0.0));	23,174✔
252
253	if (not use_initial_guess)	23,174✔
254	OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_FALSE));	20,682✔
255	else
256	OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_TRUE));	2,492✔
257
258	OpenSnPETScCall(KSPSetTolerances(	23,174✔
259	ksp_, options.residual_tolerance, options.residual_tolerance, 1.0e50, options.max_iters));
260
261	if (options.perform_symmetry_check)	23,174✔
262	{
263	PetscBool symmetry = PETSC_FALSE;	14✔
264	OpenSnPETScCall(MatIsSymmetric(A_, 1.0e-6, &symmetry));	14✔
265	if (symmetry == PETSC_FALSE)	14✔
266	throw std::logic_error(fname + ":Symmetry check failed");	×
267	}
268
269	if (options.verbose)	23,174✔
270	{
271	OpenSnPETScCall(KSPMonitorSet(ksp_, &KSPMonitorRelativeToRHS, nullptr, nullptr));	142✔
272	LogDiffusionSolveStart(name_, rhs_);	142✔
273	}
274
275	if (use_initial_guess)	23,174✔
276	{
277	PetscScalar* x_raw = nullptr;	2,492✔
278	OpenSnPETScCall(VecGetArray(x_, &x_raw));	2,492✔
279	size_t k = 0;	2,492✔
280	for (const auto& value : solution)	10,255,292✔
281	x_raw[k++] = value;	10,252,800✔
282	OpenSnPETScCall(VecRestoreArray(x_, &x_raw));	2,492✔
283	}
284
285	// Solve
286	OpenSnPETScCall(KSPSolve(ksp_, rhs_, x_));	23,174✔
287
288	// Print convergence info
289	if (options.verbose)	23,174✔
290	LogDiffusionSolveFinal(name_, ksp_, x_);	142✔
291
292	// Transfer petsc solution to vector
293	if (requires_ghosts_)	23,174✔
294	{
295	CommunicateGhostEntries(x_);	130✔
296	sdm_.LocalizePETScVectorWithGhosts(x_, solution, uk_man_);	130✔
297	}
298	else
299	sdm_.LocalizePETScVector(x_, solution, uk_man_);	23,044✔
300	}	23,174✔
301
302	void
303	DiffusionSolver::Solve(Vec petsc_solution, bool use_initial_guess)	×
304	{
305	const std::string fname = "acceleration::DiffusionMIPSolver::Solve";	×
306	OpenSnPETScCall(VecSet(x_, 0.0));	×
307
308	if (not use_initial_guess)	×
309	OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_FALSE));	×
310	else
311	OpenSnPETScCall(KSPSetInitialGuessNonzero(ksp_, PETSC_TRUE));	×
312
313	OpenSnPETScCall(KSPSetTolerances(	×
314	ksp_, options.residual_tolerance, options.residual_tolerance, 1.0e50, options.max_iters));
315
316	if (options.perform_symmetry_check)	×
317	{
318	PetscBool symmetry = PETSC_FALSE;	×
319	OpenSnPETScCall(MatIsSymmetric(A_, 1.0e-6, &symmetry));	×
320	if (symmetry == PETSC_FALSE)	×
321	throw std::logic_error(fname + ":Symmetry check failed");	×
322	}
323
324	if (options.verbose)	×
325	{
326	OpenSnPETScCall(KSPMonitorSet(ksp_, &KSPMonitorRelativeToRHS, nullptr, nullptr));	×
327	LogDiffusionSolveStart(name_, rhs_);	×
328	}
329
330	if (use_initial_guess)	×
331	{
332	OpenSnPETScCall(VecCopy(petsc_solution, x_));	×
333	}
334
335	// Solve
336	OpenSnPETScCall(KSPSolve(ksp_, rhs_, x_));	×
337
338	// Print convergence info
339	if (options.verbose)	×
340	LogDiffusionSolveFinal(name_, ksp_, x_);	×
341
342	// Transfer petsc solution to vector
343	OpenSnPETScCall(VecCopy(x_, petsc_solution));	×
344	}	×
345
346	} // namespace opensn

Open-Sn / opensn / 27259186575

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