18928565313

Committed 22 Oct 2025 07:55PM UTC coverage: 74.771%. Remained the same

Build # 18928565313

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #807 from wdhawkins/clang-tidy-init-variables

Fixing clang-tidy init-variables warnings

Run Details

132 of 179 new or added lines in 52 files covered. (73.74%)

184 existing lines in 3 files now uncovered.

18203 of 24345 relevant lines covered (74.77%)

53868061.69 hits per line

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

55.71

/modules/diffusion/diffusion_pwlc_solver.cc

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

#include "modules/diffusion/diffusion_pwlc_solver.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/acceleration/acceleration.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_structs.h"
#include "framework/math/spatial_discretization/finite_element/unit_cell_matrices.h"
#include "framework/math/spatial_discretization/spatial_discretization.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/logging/log.h"
#include "framework/utils/timer.h"
#include "framework/runtime.h"

namespace opensn
{

DiffusionPWLCSolver::DiffusionPWLCSolver(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 verbose)
  : DiffusionSolver(std::move(name),
                    sdm,
                    uk_man,
                    std::move(bcs),
                    std::move(map_mat_id_2_xs),
                    unit_cell_matrices,
                    suppress_bcs,
                    true,
                    verbose)
{
  if (sdm_.GetType() != SpatialDiscretizationType::PIECEWISE_LINEAR_CONTINUOUS)
    throw std::logic_error("acceleration::DiffusionPWLCSolver can only be used with PWLC.");
}

void
DiffusionPWLCSolver::AssembleAand_b(const std::vector<double>& q_vector)
{
  const size_t num_local_dofs = sdm_.GetNumLocalAndGhostDOFs(uk_man_);
  OpenSnInvalidArgumentIf(q_vector.size() != num_local_dofs,
                          std::string("q_vector size mismatch. ") +
                            std::to_string(q_vector.size()) + " vs " +
                            std::to_string(num_local_dofs));

  const std::string fname = "acceleration::DiffusionMIPSolver::"
                            "AssembleAand_b";
  if (A_ == nullptr or rhs_ == nullptr or ksp_ == nullptr)
    throw std::logic_error(fname + ": Some or all PETSc elements are null. "
                                   "Check that Initialize has been called.");
  if (options.verbose)
    log.Log() << program_timer.GetTimeString() << " Starting assembly";

  const size_t num_groups = uk_man_.unknowns.front().num_components;

  VecSet(rhs_, 0.0);
  for (const auto& cell : grid_->local_cells)
  {
    const size_t num_faces = cell.faces.size();
    const auto& cell_mapping = sdm_.GetCellMapping(cell);
    const auto num_nodes = cell_mapping.GetNumNodes();
    const auto cc_nodes = cell_mapping.GetNodeLocations();
    const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id];

    const auto& intV_gradshapeI_gradshapeJ = unit_cell_matrices.intV_gradshapeI_gradshapeJ;
    const auto& intV_shapeI_shapeJ = unit_cell_matrices.intV_shapeI_shapeJ;

    const auto& xs = mat_id_2_xs_map_.at(cell.block_id);

    // Mark dirichlet nodes
    std::vector<std::pair<bool, double>> node_is_dirichlet(num_nodes, {false, 0.0});
    for (size_t f = 0; f < num_faces; ++f)
    {
      const auto& face = cell.faces[f];
      if (not face.has_neighbor and not suppress_bcs_)
      {
        BoundaryCondition bc;
        if (bcs_.count(face.neighbor_id) > 0)
          bc = bcs_.at(face.neighbor_id);

        if (bc.type != BCType::DIRICHLET)
          continue;

        const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);
        for (size_t fi = 0; fi < num_face_nodes; ++fi)
          node_is_dirichlet[cell_mapping.MapFaceNode(f, fi)] = {true, bc.values[0]};
      }
    }

    DenseMatrix<double> cell_A(num_nodes, num_nodes);
    Vector<double> cell_rhs(num_nodes);
    Vector<PetscInt> cell_idxs(num_nodes);
    for (size_t g = 0; g < num_groups; ++g)
    {
      cell_A.Set(0.);
      cell_rhs.Set(0.);
      for (size_t i = 0; i < num_nodes; ++i)
        cell_idxs(i) = static_cast<PetscInt>(sdm_.MapDOF(cell, i, uk_man_, 0, g));

      // Get coefficient and nodal src
      const double Dg = xs.Dg[g];
      const double sigr_g = xs.sigR[g];

      std::vector<double> qg(num_nodes, 0.0);
      for (size_t j = 0; j < num_nodes; ++j)
        qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];

      // Assemble continuous terms
      for (size_t i = 0; i < num_nodes; ++i)
      {
        if (node_is_dirichlet[i].first)
          continue;
        double entry_rhs_i = 0.0;
        for (size_t j = 0; j < num_nodes; ++j)
        {

          const double entry_aij =
            Dg * intV_gradshapeI_gradshapeJ(i, j) + sigr_g * intV_shapeI_shapeJ(i, j);

          if (not node_is_dirichlet[j].first)
            cell_A(i, j) += entry_aij;
          else
          {
            const double bcvalue = node_is_dirichlet[j].second;
            cell_rhs(i) -= entry_aij * bcvalue;
          }

          entry_rhs_i += intV_shapeI_shapeJ(i, j) * qg[j];
        } // for j

        cell_rhs(i) = entry_rhs_i;
      } // for i

      // Assemble face terms
      for (size_t f = 0; f < num_faces; ++f)
      {
        const auto& face = cell.faces[f];
        const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);

        const auto& intS_shapeI_shapeJ = unit_cell_matrices.intS_shapeI_shapeJ[f];
        const auto& intS_shapeI = unit_cell_matrices.intS_shapeI[f];

        if (not face.has_neighbor and not suppress_bcs_)
        {
          BoundaryCondition bc;
          if (bcs_.count(face.neighbor_id) > 0)
            bc = bcs_.at(face.neighbor_id);

          if (bc.type == BCType::DIRICHLET)
          {
            const double bc_value = bc.values[0];

            for (size_t fi = 0; fi < num_face_nodes; ++fi)
            {
              const int i = cell_mapping.MapFaceNode(f, fi);

              // MatSetValue(A_, imap, imap, intV_shapeI[i], ADD_VALUES);
              // VecSetValue(rhs_, imap, bc_value * intV_shapeI[i], ADD_VALUES);
              cell_A(i, i) += 1.0;
              cell_rhs(i) += bc_value;
            } // for fi

          } // Dirichlet BC
          else if (bc.type == BCType::ROBIN)
          {
            const double aval = bc.values[0];
            const double bval = bc.values[1];
            const double fval = bc.values[2];

            if (std::fabs(bval) < 1.0e-12)
              continue; // a and f assumed zero

            for (size_t fi = 0; fi < num_face_nodes; ++fi)
            {
              const int i = cell_mapping.MapFaceNode(f, fi);

              if (std::fabs(aval) >= 1.0e-12)
              {
                for (size_t fj = 0; fj < num_face_nodes; ++fj)
                {
                  const int j = cell_mapping.MapFaceNode(f, fj);

                  const double aij = (aval / bval) * intS_shapeI_shapeJ(i, j);

                  cell_A(i, j) += aij;
                } // for fj
              } // if a nonzero

              if (std::fabs(fval) >= 1.0e-12)
              {
                const double rhs_val = (fval / bval) * intS_shapeI(i);

                cell_rhs(i) += rhs_val;
              } // if f nonzero
            } // for fi
          } // Robin BC
        } // boundary face
      } // for face

      auto nn = static_cast<PetscInt>(num_nodes);
      MatSetValues(A_, nn, cell_idxs.data(), nn, cell_idxs.data(), cell_A.data(), ADD_VALUES);
      VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);
    } // for g
  } // for cell

  MatAssemblyBegin(A_, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(A_, MAT_FINAL_ASSEMBLY);
  VecAssemblyBegin(rhs_);
  VecAssemblyEnd(rhs_);

  if (options.verbose)
  {
    MatInfo info;
    MatGetInfo(A_, MAT_GLOBAL_SUM, &info);

    log.Log() << "Number of mallocs used = " << info.mallocs
              << "\nNumber of non-zeros allocated = " << info.nz_allocated
              << "\nNumber of non-zeros used = " << info.nz_used
              << "\nNumber of unneeded non-zeros = " << info.nz_unneeded;
  }

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

  KSPSetOperators(ksp_, A_, A_);

  if (options.verbose)
    log.Log() << program_timer.GetTimeString() << " Assembly completed";

  PC pc = nullptr;
  KSPGetPC(ksp_, &pc);
  PCSetUp(pc);

  KSPSetUp(ksp_);
}

void
DiffusionPWLCSolver::Assemble_b(const std::vector<double>& q_vector)
{
  const size_t num_local_dofs = sdm_.GetNumLocalAndGhostDOFs(uk_man_);
  OpenSnInvalidArgumentIf(q_vector.size() != num_local_dofs,
                          std::string("q_vector size mismatch. ") +
                            std::to_string(q_vector.size()) + " vs " +
                            std::to_string(num_local_dofs));
  const std::string fname = "acceleration::DiffusionMIPSolver::"
                            "Assemble_b";
  if (A_ == nullptr or rhs_ == nullptr or ksp_ == nullptr)
    throw std::logic_error(fname + ": Some or all PETSc elements are null. "
                                   "Check that Initialize has been called.");
  if (options.verbose)
    log.Log() << program_timer.GetTimeString() << " Starting assembly";

  const size_t num_groups = uk_man_.unknowns.front().num_components;

  VecSet(rhs_, 0.0);
  for (const auto& cell : grid_->local_cells)
  {
    const size_t num_faces = cell.faces.size();
    const auto& cell_mapping = sdm_.GetCellMapping(cell);
    const auto num_nodes = cell_mapping.GetNumNodes();
    const auto cc_nodes = cell_mapping.GetNodeLocations();
    const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id];

    const auto& intV_gradshapeI_gradshapeJ = unit_cell_matrices.intV_gradshapeI_gradshapeJ;
    const auto& intV_shapeI_shapeJ = unit_cell_matrices.intV_shapeI_shapeJ;
    const auto& intV_shapeI = unit_cell_matrices.intV_shapeI;

    const auto& xs = mat_id_2_xs_map_.at(cell.block_id);

    // Mark dirichlet nodes
    std::vector<std::pair<bool, double>> node_is_dirichlet(num_nodes, {false, 0.0});
    for (size_t f = 0; f < num_faces; ++f)
    {
      const auto& face = cell.faces[f];
      if (not face.has_neighbor and suppress_bcs_)
      {
        BoundaryCondition bc;
        if (bcs_.count(face.neighbor_id) > 0)
          bc = bcs_.at(face.neighbor_id);

        if (bc.type != BCType::DIRICHLET)
          continue;

        const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);
        for (size_t fi = 0; fi < num_face_nodes; ++fi)
          node_is_dirichlet[cell_mapping.MapFaceNode(f, fi)] = {true, bc.values[0]};
      }
    }

    Vector<double> cell_rhs(num_nodes);
    Vector<PetscInt> cell_idxs(num_nodes);

    for (size_t g = 0; g < num_groups; ++g)
    {
      cell_rhs.Set(0.);
      for (size_t i = 0; i < num_nodes; ++i)
        cell_idxs(i) = static_cast<PetscInt>(sdm_.MapDOF(cell, i, uk_man_, 0, g));

      // Get coefficient and nodal src
      std::vector<double> qg(num_nodes, 0.0);
      for (size_t j = 0; j < num_nodes; ++j)
        qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];

      // Assemble continuous terms
      const double Dg = xs.Dg[g];
      const double sigr_g = xs.sigR[g];

      for (size_t i = 0; i < num_nodes; ++i)
      {
        if (node_is_dirichlet[i].first)
          continue;
        double entry_rhs_i = 0.0;
        for (size_t j = 0; j < num_nodes; ++j)
        {
          if (node_is_dirichlet[j].first)
          {
            const double entry_aij =
              Dg * intV_gradshapeI_gradshapeJ(i, j) + sigr_g * intV_shapeI_shapeJ(i, j);

            const double bcvalue = node_is_dirichlet[j].second;
            cell_rhs(i) -= entry_aij * bcvalue;
          }

          entry_rhs_i += intV_shapeI_shapeJ(i, j) * qg[j];
        } // for j

        cell_rhs(i) += entry_rhs_i;
      } // for i

      // Assemble face terms
      for (size_t f = 0; f < num_faces; ++f)
      {
        const auto& face = cell.faces[f];
        const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);

        const auto& intS_shapeI = unit_cell_matrices.intS_shapeI[f];

        if (not face.has_neighbor and suppress_bcs_)
        {
          BoundaryCondition bc;
          if (bcs_.count(face.neighbor_id) > 0)
            bc = bcs_.at(face.neighbor_id);

          if (bc.type == BCType::DIRICHLET)
          {
            const double bc_value = bc.values[0];

            // Assembly penalty terms
            for (size_t fi = 0; fi < num_face_nodes; ++fi)
            {
              const int i = cell_mapping.MapFaceNode(f, fi);

              // VecSetValue(rhs_, imap, bc_value * intV_shapeI[i], ADD_VALUES);
              cell_rhs(i) += bc_value;
            } // for fi

          } // Dirichlet BC
          else if (bc.type == BCType::ROBIN)
          {
            const double bval = bc.values[1];
            const double fval = bc.values[2];

            if (std::fabs(bval) < 1.0e-12)
              continue; // a and f assumed zero

            for (size_t fi = 0; fi < num_face_nodes; ++fi)
            {
              const int i = cell_mapping.MapFaceNode(f, fi);

              if (std::fabs(fval) >= 1.0e-12)
              {
                const double rhs_val = (fval / bval) * intS_shapeI(i);

                cell_rhs(i) += rhs_val;
              } // if f nonzero
            } // for fi
          } // Robin BC
        } // boundary face
      } // for face
      auto nn = static_cast<PetscInt>(num_nodes);
      VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);
    } // for g
  } // for cell

  VecAssemblyBegin(rhs_);
  VecAssemblyEnd(rhs_);

  if (options.verbose)
    log.Log() << program_timer.GetTimeString() << " Assembly completed";
}

void
DiffusionPWLCSolver::Assemble_b(Vec petsc_q_vector)
{
  const std::string fname = "acceleration::DiffusionMIPSolver::"
                            "Assemble_b";
  if (A_ == nullptr or rhs_ == nullptr or ksp_ == nullptr)
    throw std::logic_error(fname + ": Some or all PETSc elements are null. "
                                   "Check that Initialize has been called.");
  if (options.verbose)
    log.Log() << program_timer.GetTimeString() << " Starting assembly";

  const size_t num_groups = uk_man_.unknowns.front().num_components;

  const double* q_vector = nullptr;
  VecGetArrayRead(petsc_q_vector, &q_vector);

  VecSet(rhs_, 0.0);
  for (const auto& cell : grid_->local_cells)
  {
    const size_t num_faces = cell.faces.size();
    const auto& cell_mapping = sdm_.GetCellMapping(cell);
    const auto num_nodes = cell_mapping.GetNumNodes();
    const auto cc_nodes = cell_mapping.GetNodeLocations();
    const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id];

    const auto& intV_shapeI_shapeJ = unit_cell_matrices.intV_shapeI_shapeJ;
    const auto& intV_shapeI = unit_cell_matrices.intV_shapeI;

    // Mark dirichlet nodes
    std::vector<bool> node_is_dirichlet(num_nodes, false);
    for (size_t f = 0; f < num_faces; ++f)
    {
      const auto& face = cell.faces[f];
      if (not face.has_neighbor and not suppress_bcs_)
      {
        BoundaryCondition bc;
        if (bcs_.count(face.neighbor_id) > 0)
          bc = bcs_.at(face.neighbor_id);

        if (bc.type != BCType::DIRICHLET)
          continue;

        const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);
        for (size_t fi = 0; fi < num_face_nodes; ++fi)
          node_is_dirichlet[cell_mapping.MapFaceNode(f, fi)] = true;
      }
    }

    Vector<double> cell_rhs(num_nodes);
    Vector<PetscInt> cell_idxs(num_nodes);

    for (size_t g = 0; g < num_groups; ++g)
    {
      cell_rhs.Set(0.);
      for (size_t i = 0; i < num_nodes; ++i)
        cell_idxs(i) = static_cast<PetscInt>(sdm_.MapDOF(cell, i, uk_man_, 0, g));

      // Get coefficient and nodal src
      std::vector<double> qg(num_nodes, 0.0);
      for (size_t j = 0; j < num_nodes; ++j)
        qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];

      // Assemble continuous terms
      for (size_t i = 0; i < num_nodes; ++i)
      {
        if (node_is_dirichlet[i])
          continue;
        double entry_rhs_i = 0.0;
        for (size_t j = 0; j < num_nodes; ++j)
          entry_rhs_i += intV_shapeI_shapeJ(i, j) * qg[j];

        cell_rhs(i) += entry_rhs_i;
      } // for i

      // Assemble face terms
      for (size_t f = 0; f < num_faces; ++f)
      {
        const auto& face = cell.faces[f];
        const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);

        const auto& intS_shapeI = unit_cell_matrices.intS_shapeI[f];

        if (not face.has_neighbor and not suppress_bcs_)
        {
          BoundaryCondition bc;
          if (bcs_.count(face.neighbor_id) > 0)
            bc = bcs_.at(face.neighbor_id);

          if (bc.type == BCType::DIRICHLET)
          {
            const double bc_value = bc.values[0];

            // Assembly penalty terms
            for (size_t fi = 0; fi < num_face_nodes; ++fi)
            {
              const int i = cell_mapping.MapFaceNode(f, fi);

              cell_rhs(i) += bc_value * intV_shapeI(i);
            } // for fi

          } // Dirichlet BC
          else if (bc.type == BCType::ROBIN)
          {
            const double bval = bc.values[1];
            const double fval = bc.values[2];

            if (std::fabs(bval) < 1.0e-12)
              continue; // a and f assumed zero

            for (size_t fi = 0; fi < num_face_nodes; ++fi)
            {
              const int i = cell_mapping.MapFaceNode(f, fi);

              if (std::fabs(fval) >= 1.0e-12)
              {
                const double rhs_val = (fval / bval) * intS_shapeI(i);

                cell_rhs(i) += rhs_val;
              } // if f nonzero
            } // for fi
          } // Robin BC
        } // boundary face
      } // for face
      auto nn = static_cast<PetscInt>(num_nodes);
      VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);
    } // for g
  } // for cell

  VecRestoreArrayRead(petsc_q_vector, &q_vector);

  VecAssemblyBegin(rhs_);
  VecAssemblyEnd(rhs_);

  if (options.verbose)
    log.Log() << program_timer.GetTimeString() << " Assembly completed";
}

} // 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_pwlc_solver.h"
5	#include "modules/linear_boltzmann_solvers/lbs_problem/acceleration/acceleration.h"
6	#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_structs.h"
7	#include "framework/math/spatial_discretization/finite_element/unit_cell_matrices.h"
8	#include "framework/math/spatial_discretization/spatial_discretization.h"
9	#include "framework/mesh/mesh_continuum/mesh_continuum.h"
10	#include "framework/logging/log.h"
11	#include "framework/utils/timer.h"
12	#include "framework/runtime.h"
13
14	namespace opensn
15	{
16
17	DiffusionPWLCSolver::DiffusionPWLCSolver(std::string name,	1✔
18	const opensn::SpatialDiscretization& sdm,
19	const UnknownManager& uk_man,
20	std::map<uint64_t, BoundaryCondition> bcs,
21	MatID2XSMap map_mat_id_2_xs,
22	const std::vector<UnitCellMatrices>& unit_cell_matrices,
23	const bool suppress_bcs,
24	const bool verbose)	1✔
25	: DiffusionSolver(std::move(name),
26	sdm,
27	uk_man,
28	std::move(bcs),
29	std::move(map_mat_id_2_xs),
30	unit_cell_matrices,
31	suppress_bcs,
32	true,
33	verbose)	2✔
34	{
35	if (sdm_.GetType() != SpatialDiscretizationType::PIECEWISE_LINEAR_CONTINUOUS)	1✔
36	throw std::logic_error("acceleration::DiffusionPWLCSolver can only be used with PWLC.");	×
37	}	1✔
38
39	void
40	DiffusionPWLCSolver::AssembleAand_b(const std::vector<double>& q_vector)	1✔
41	{
42	const size_t num_local_dofs = sdm_.GetNumLocalAndGhostDOFs(uk_man_);	1✔
43	OpenSnInvalidArgumentIf(q_vector.size() != num_local_dofs,	1✔
44	std::string("q_vector size mismatch. ") +
45	std::to_string(q_vector.size()) + " vs " +
46	std::to_string(num_local_dofs));
47
48	const std::string fname = "acceleration::DiffusionMIPSolver::"	1✔
49	"AssembleAand_b";	1✔
50	if (A_ == nullptr or rhs_ == nullptr or ksp_ == nullptr)	1✔
51	throw std::logic_error(fname + ": Some or all PETSc elements are null. "	×
52	"Check that Initialize has been called.");	×
53	if (options.verbose)	1✔
54	log.Log() << program_timer.GetTimeString() << " Starting assembly";	3✔
55
56	const size_t num_groups = uk_man_.unknowns.front().num_components;	1✔
57
58	VecSet(rhs_, 0.0);	1✔
59	for (const auto& cell : grid_->local_cells)	101✔
60	{
61	const size_t num_faces = cell.faces.size();	100✔
62	const auto& cell_mapping = sdm_.GetCellMapping(cell);	100✔
63	const auto num_nodes = cell_mapping.GetNumNodes();	100✔
64	const auto cc_nodes = cell_mapping.GetNodeLocations();	100✔
65	const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id];	100✔
66
67	const auto& intV_gradshapeI_gradshapeJ = unit_cell_matrices.intV_gradshapeI_gradshapeJ;	100✔
68	const auto& intV_shapeI_shapeJ = unit_cell_matrices.intV_shapeI_shapeJ;	100✔
69
70	const auto& xs = mat_id_2_xs_map_.at(cell.block_id);	100✔
71
72	// Mark dirichlet nodes
73	std::vector<std::pair<bool, double>> node_is_dirichlet(num_nodes, {false, 0.0});	100✔
74	for (size_t f = 0; f < num_faces; ++f)	500✔
75	{
76	const auto& face = cell.faces[f];	400✔
77	if (not face.has_neighbor and not suppress_bcs_)	400✔
78	{
79	BoundaryCondition bc;	40✔
80	if (bcs_.count(face.neighbor_id) > 0)	40✔
81	bc = bcs_.at(face.neighbor_id);	40✔
82
83	if (bc.type != BCType::DIRICHLET)	40✔
84	continue;	40✔
85
86	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	×
87	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
88	node_is_dirichlet[cell_mapping.MapFaceNode(f, fi)] = {true, bc.values[0]};	×
89	}
90	}
91
92	DenseMatrix<double> cell_A(num_nodes, num_nodes);	100✔
93	Vector<double> cell_rhs(num_nodes);	100✔
94	Vector<PetscInt> cell_idxs(num_nodes);	100✔
95	for (size_t g = 0; g < num_groups; ++g)	200✔
96	{
97	cell_A.Set(0.);	100✔
98	cell_rhs.Set(0.);	100✔
99	for (size_t i = 0; i < num_nodes; ++i)	500✔
100	cell_idxs(i) = static_cast<PetscInt>(sdm_.MapDOF(cell, i, uk_man_, 0, g));	400✔
101
102	// Get coefficient and nodal src
103	const double Dg = xs.Dg[g];	100✔
104	const double sigr_g = xs.sigR[g];	100✔
105
106	std::vector<double> qg(num_nodes, 0.0);	100✔
107	for (size_t j = 0; j < num_nodes; ++j)	500✔
108	qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];	400✔
109
110	// Assemble continuous terms
111	for (size_t i = 0; i < num_nodes; ++i)	500✔
112	{
113	if (node_is_dirichlet[i].first)	400✔
114	continue;	×
115	double entry_rhs_i = 0.0;
116	for (size_t j = 0; j < num_nodes; ++j)	2,000✔
117	{
118
119	const double entry_aij =	1,600✔
120	Dg * intV_gradshapeI_gradshapeJ(i, j) + sigr_g * intV_shapeI_shapeJ(i, j);	1,600✔
121
122	if (not node_is_dirichlet[j].first)	1,600✔
123	cell_A(i, j) += entry_aij;	1,600✔
124	else
125	{
126	const double bcvalue = node_is_dirichlet[j].second;	×
127	cell_rhs(i) -= entry_aij * bcvalue;	×
128	}
129
130	entry_rhs_i += intV_shapeI_shapeJ(i, j) * qg[j];	1,600✔
131	} // for j
132
133	cell_rhs(i) = entry_rhs_i;	400✔
134	} // for i
135
136	// Assemble face terms
137	for (size_t f = 0; f < num_faces; ++f)	500✔
138	{
139	const auto& face = cell.faces[f];	400✔
140	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	400✔
141
142	const auto& intS_shapeI_shapeJ = unit_cell_matrices.intS_shapeI_shapeJ[f];	400✔
143	const auto& intS_shapeI = unit_cell_matrices.intS_shapeI[f];	400✔
144
145	if (not face.has_neighbor and not suppress_bcs_)	400✔
146	{
147	BoundaryCondition bc;	40✔
148	if (bcs_.count(face.neighbor_id) > 0)	40✔
149	bc = bcs_.at(face.neighbor_id);	40✔
150
151	if (bc.type == BCType::DIRICHLET)	40✔
152	{
153	const double bc_value = bc.values[0];
154
155	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
156	{
157	const int i = cell_mapping.MapFaceNode(f, fi);	×
158
159	// MatSetValue(A_, imap, imap, intV_shapeI[i], ADD_VALUES);
160	// VecSetValue(rhs_, imap, bc_value * intV_shapeI[i], ADD_VALUES);
161	cell_A(i, i) += 1.0;	×
162	cell_rhs(i) += bc_value;	×
163	} // for fi
164
165	} // Dirichlet BC
166	else if (bc.type == BCType::ROBIN)	40✔
167	{
168	const double aval = bc.values[0];	40✔
169	const double bval = bc.values[1];	40✔
170	const double fval = bc.values[2];	40✔
171
172	if (std::fabs(bval) < 1.0e-12)	40✔
173	continue; // a and f assumed zero	×
174
175	for (size_t fi = 0; fi < num_face_nodes; ++fi)	120✔
176	{
177	const int i = cell_mapping.MapFaceNode(f, fi);	80✔
178
179	if (std::fabs(aval) >= 1.0e-12)	80✔
180	{
181	for (size_t fj = 0; fj < num_face_nodes; ++fj)	240✔
182	{
183	const int j = cell_mapping.MapFaceNode(f, fj);	160✔
184
185	const double aij = (aval / bval) * intS_shapeI_shapeJ(i, j);	160✔
186
187	cell_A(i, j) += aij;	160✔
188	} // for fj
189	} // if a nonzero
190
191	if (std::fabs(fval) >= 1.0e-12)	80✔
192	{
193	const double rhs_val = (fval / bval) * intS_shapeI(i);	×
194
195	cell_rhs(i) += rhs_val;	×
196	} // if f nonzero
197	} // for fi
198	} // Robin BC
199	} // boundary face
200	} // for face
201
202	auto nn = static_cast<PetscInt>(num_nodes);	100✔
203	MatSetValues(A_, nn, cell_idxs.data(), nn, cell_idxs.data(), cell_A.data(), ADD_VALUES);	100✔
204	VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);	100✔
205	} // for g	100✔
206	} // for cell	300✔
207
208	MatAssemblyBegin(A_, MAT_FINAL_ASSEMBLY);	1✔
209	MatAssemblyEnd(A_, MAT_FINAL_ASSEMBLY);	1✔
210	VecAssemblyBegin(rhs_);	1✔
211	VecAssemblyEnd(rhs_);	1✔
212
213	if (options.verbose)	1✔
214	{
215	MatInfo info;	1✔
216	MatGetInfo(A_, MAT_GLOBAL_SUM, &info);	1✔
217
218	log.Log() << "Number of mallocs used = " << info.mallocs	3✔
219	<< "\nNumber of non-zeros allocated = " << info.nz_allocated	1✔
220	<< "\nNumber of non-zeros used = " << info.nz_used	1✔
221	<< "\nNumber of unneeded non-zeros = " << info.nz_unneeded;	2✔
222	}
223
224	if (options.perform_symmetry_check)	1✔
225	{
226	PetscBool symmetry = PETSC_FALSE;	1✔
227	MatIsSymmetric(A_, 1.0e-6, &symmetry);	1✔
228	if (symmetry == PETSC_FALSE)	1✔
229	throw std::logic_error(fname + ":Symmetry check failed");	×
230	}
231
232	KSPSetOperators(ksp_, A_, A_);	1✔
233
234	if (options.verbose)	1✔
235	log.Log() << program_timer.GetTimeString() << " Assembly completed";	3✔
236
237	PC pc = nullptr;	1✔
238	KSPGetPC(ksp_, &pc);	1✔
239	PCSetUp(pc);	1✔
240
241	KSPSetUp(ksp_);	1✔
242	}	1✔
243
244	void
245	DiffusionPWLCSolver::Assemble_b(const std::vector<double>& q_vector)	1✔
246	{
247	const size_t num_local_dofs = sdm_.GetNumLocalAndGhostDOFs(uk_man_);	1✔
248	OpenSnInvalidArgumentIf(q_vector.size() != num_local_dofs,	1✔
249	std::string("q_vector size mismatch. ") +
250	std::to_string(q_vector.size()) + " vs " +
251	std::to_string(num_local_dofs));
252	const std::string fname = "acceleration::DiffusionMIPSolver::"	1✔
253	"Assemble_b";	1✔
254	if (A_ == nullptr or rhs_ == nullptr or ksp_ == nullptr)	1✔
255	throw std::logic_error(fname + ": Some or all PETSc elements are null. "	×
256	"Check that Initialize has been called.");	×
257	if (options.verbose)	1✔
258	log.Log() << program_timer.GetTimeString() << " Starting assembly";	3✔
259
260	const size_t num_groups = uk_man_.unknowns.front().num_components;	1✔
261
262	VecSet(rhs_, 0.0);	1✔
263	for (const auto& cell : grid_->local_cells)	101✔
264	{
265	const size_t num_faces = cell.faces.size();	100✔
266	const auto& cell_mapping = sdm_.GetCellMapping(cell);	100✔
267	const auto num_nodes = cell_mapping.GetNumNodes();	100✔
268	const auto cc_nodes = cell_mapping.GetNodeLocations();	100✔
269	const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id];	100✔
270
271	const auto& intV_gradshapeI_gradshapeJ = unit_cell_matrices.intV_gradshapeI_gradshapeJ;	100✔
272	const auto& intV_shapeI_shapeJ = unit_cell_matrices.intV_shapeI_shapeJ;	100✔
273	const auto& intV_shapeI = unit_cell_matrices.intV_shapeI;	100✔
274
275	const auto& xs = mat_id_2_xs_map_.at(cell.block_id);	100✔
276
277	// Mark dirichlet nodes
278	std::vector<std::pair<bool, double>> node_is_dirichlet(num_nodes, {false, 0.0});	100✔
279	for (size_t f = 0; f < num_faces; ++f)	500✔
280	{
281	const auto& face = cell.faces[f];	400✔
282	if (not face.has_neighbor and suppress_bcs_)	400✔
283	{
284	BoundaryCondition bc;	×
285	if (bcs_.count(face.neighbor_id) > 0)	×
286	bc = bcs_.at(face.neighbor_id);	×
287
288	if (bc.type != BCType::DIRICHLET)	×
289	continue;	×
290
291	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	×
292	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
293	node_is_dirichlet[cell_mapping.MapFaceNode(f, fi)] = {true, bc.values[0]};	×
294	}
295	}
296
297	Vector<double> cell_rhs(num_nodes);	100✔
298	Vector<PetscInt> cell_idxs(num_nodes);	100✔
299
300	for (size_t g = 0; g < num_groups; ++g)	200✔
301	{
302	cell_rhs.Set(0.);	100✔
303	for (size_t i = 0; i < num_nodes; ++i)	500✔
304	cell_idxs(i) = static_cast<PetscInt>(sdm_.MapDOF(cell, i, uk_man_, 0, g));	400✔
305
306	// Get coefficient and nodal src
307	std::vector<double> qg(num_nodes, 0.0);	100✔
308	for (size_t j = 0; j < num_nodes; ++j)	500✔
309	qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];	400✔
310
311	// Assemble continuous terms
312	const double Dg = xs.Dg[g];	100✔
313	const double sigr_g = xs.sigR[g];	100✔
314
315	for (size_t i = 0; i < num_nodes; ++i)	500✔
316	{
317	if (node_is_dirichlet[i].first)	400✔
318	continue;	×
319	double entry_rhs_i = 0.0;
320	for (size_t j = 0; j < num_nodes; ++j)	2,000✔
321	{
322	if (node_is_dirichlet[j].first)	1,600✔
323	{
324	const double entry_aij =	×
325	Dg * intV_gradshapeI_gradshapeJ(i, j) + sigr_g * intV_shapeI_shapeJ(i, j);	×
326
327	const double bcvalue = node_is_dirichlet[j].second;	×
328	cell_rhs(i) -= entry_aij * bcvalue;	×
329	}
330
331	entry_rhs_i += intV_shapeI_shapeJ(i, j) * qg[j];	1,600✔
332	} // for j
333
334	cell_rhs(i) += entry_rhs_i;	400✔
335	} // for i
336
337	// Assemble face terms
338	for (size_t f = 0; f < num_faces; ++f)	500✔
339	{
340	const auto& face = cell.faces[f];	400✔
341	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	400✔
342
343	const auto& intS_shapeI = unit_cell_matrices.intS_shapeI[f];	400✔
344
345	if (not face.has_neighbor and suppress_bcs_)	400✔
346	{
347	BoundaryCondition bc;	×
348	if (bcs_.count(face.neighbor_id) > 0)	×
349	bc = bcs_.at(face.neighbor_id);	×
350
351	if (bc.type == BCType::DIRICHLET)	×
352	{
353	const double bc_value = bc.values[0];
354
355	// Assembly penalty terms
356	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
357	{
358	const int i = cell_mapping.MapFaceNode(f, fi);	×
359
360	// VecSetValue(rhs_, imap, bc_value * intV_shapeI[i], ADD_VALUES);
361	cell_rhs(i) += bc_value;	×
362	} // for fi
363
364	} // Dirichlet BC
365	else if (bc.type == BCType::ROBIN)	×
366	{
367	const double bval = bc.values[1];	×
368	const double fval = bc.values[2];	×
369
370	if (std::fabs(bval) < 1.0e-12)	×
371	continue; // a and f assumed zero	×
372
373	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
374	{
375	const int i = cell_mapping.MapFaceNode(f, fi);	×
376
377	if (std::fabs(fval) >= 1.0e-12)	×
378	{
379	const double rhs_val = (fval / bval) * intS_shapeI(i);	×
380
381	cell_rhs(i) += rhs_val;	×
382	} // if f nonzero
383	} // for fi
384	} // Robin BC
385	} // boundary face
386	} // for face
387	auto nn = static_cast<PetscInt>(num_nodes);	100✔
388	VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);	100✔
389	} // for g	100✔
390	} // for cell	200✔
391
392	VecAssemblyBegin(rhs_);	1✔
393	VecAssemblyEnd(rhs_);	1✔
394
395	if (options.verbose)	1✔
396	log.Log() << program_timer.GetTimeString() << " Assembly completed";	3✔
397	}	1✔
398
399	void
400	DiffusionPWLCSolver::Assemble_b(Vec petsc_q_vector)	×
401	{
402	const std::string fname = "acceleration::DiffusionMIPSolver::"	×
403	"Assemble_b";	×
404	if (A_ == nullptr or rhs_ == nullptr or ksp_ == nullptr)	×
405	throw std::logic_error(fname + ": Some or all PETSc elements are null. "	×
406	"Check that Initialize has been called.");	×
407	if (options.verbose)	×
408	log.Log() << program_timer.GetTimeString() << " Starting assembly";	×
409
410	const size_t num_groups = uk_man_.unknowns.front().num_components;	×
411
NEW 412	const double* q_vector = nullptr;	×
413	VecGetArrayRead(petsc_q_vector, &q_vector);	×
414
415	VecSet(rhs_, 0.0);	×
416	for (const auto& cell : grid_->local_cells)	×
417	{
418	const size_t num_faces = cell.faces.size();	×
419	const auto& cell_mapping = sdm_.GetCellMapping(cell);	×
420	const auto num_nodes = cell_mapping.GetNumNodes();	×
421	const auto cc_nodes = cell_mapping.GetNodeLocations();	×
422	const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id];	×
423
424	const auto& intV_shapeI_shapeJ = unit_cell_matrices.intV_shapeI_shapeJ;	×
425	const auto& intV_shapeI = unit_cell_matrices.intV_shapeI;	×
426
427	// Mark dirichlet nodes
428	std::vector<bool> node_is_dirichlet(num_nodes, false);	×
429	for (size_t f = 0; f < num_faces; ++f)	×
430	{
431	const auto& face = cell.faces[f];	×
432	if (not face.has_neighbor and not suppress_bcs_)	×
433	{
434	BoundaryCondition bc;	×
435	if (bcs_.count(face.neighbor_id) > 0)	×
436	bc = bcs_.at(face.neighbor_id);	×
437
438	if (bc.type != BCType::DIRICHLET)	×
439	continue;	×
440
441	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	×
442	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
443	node_is_dirichlet[cell_mapping.MapFaceNode(f, fi)] = true;	×
444	}
445	}
446
447	Vector<double> cell_rhs(num_nodes);	×
448	Vector<PetscInt> cell_idxs(num_nodes);	×
449
450	for (size_t g = 0; g < num_groups; ++g)	×
451	{
452	cell_rhs.Set(0.);	×
453	for (size_t i = 0; i < num_nodes; ++i)	×
454	cell_idxs(i) = static_cast<PetscInt>(sdm_.MapDOF(cell, i, uk_man_, 0, g));	×
455
456	// Get coefficient and nodal src
457	std::vector<double> qg(num_nodes, 0.0);	×
458	for (size_t j = 0; j < num_nodes; ++j)	×
459	qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];	×
460
461	// Assemble continuous terms
462	for (size_t i = 0; i < num_nodes; ++i)	×
463	{
464	if (node_is_dirichlet[i])	×
465	continue;	×
466	double entry_rhs_i = 0.0;
467	for (size_t j = 0; j < num_nodes; ++j)	×
468	entry_rhs_i += intV_shapeI_shapeJ(i, j) * qg[j];	×
469
470	cell_rhs(i) += entry_rhs_i;	×
471	} // for i
472
473	// Assemble face terms
474	for (size_t f = 0; f < num_faces; ++f)	×
475	{
476	const auto& face = cell.faces[f];	×
477	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	×
478
479	const auto& intS_shapeI = unit_cell_matrices.intS_shapeI[f];	×
480
481	if (not face.has_neighbor and not suppress_bcs_)	×
482	{
483	BoundaryCondition bc;	×
484	if (bcs_.count(face.neighbor_id) > 0)	×
485	bc = bcs_.at(face.neighbor_id);	×
486
487	if (bc.type == BCType::DIRICHLET)	×
488	{
489	const double bc_value = bc.values[0];
490
491	// Assembly penalty terms
492	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
493	{
494	const int i = cell_mapping.MapFaceNode(f, fi);	×
495
496	cell_rhs(i) += bc_value * intV_shapeI(i);	×
497	} // for fi
498
499	} // Dirichlet BC
500	else if (bc.type == BCType::ROBIN)	×
501	{
502	const double bval = bc.values[1];	×
503	const double fval = bc.values[2];	×
504
505	if (std::fabs(bval) < 1.0e-12)	×
506	continue; // a and f assumed zero	×
507
508	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
509	{
510	const int i = cell_mapping.MapFaceNode(f, fi);	×
511
512	if (std::fabs(fval) >= 1.0e-12)	×
513	{
514	const double rhs_val = (fval / bval) * intS_shapeI(i);	×
515
516	cell_rhs(i) += rhs_val;	×
517	} // if f nonzero
518	} // for fi
519	} // Robin BC
520	} // boundary face
521	} // for face
522	auto nn = static_cast<PetscInt>(num_nodes);	×
523	VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);	×
524	} // for g	×
525	} // for cell	×
526
527	VecRestoreArrayRead(petsc_q_vector, &q_vector);	×
528
529	VecAssemblyBegin(rhs_);	×
530	VecAssemblyEnd(rhs_);	×
531
532	if (options.verbose)	×
533	log.Log() << program_timer.GetTimeString() << " Assembly completed";	×
534	}	×
535
536	} // namespace opensn

Open-Sn / opensn / 18928565313

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