22249939997

Committed 20 Feb 2026 06:11PM UTC coverage: 74.156% (-0.2%) from 74.335%

Build # 22249939997

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #937 from andrsd/issue/90-num-precursors

Number of precursors is `unsigned int`

Run Details

4 of 7 new or added lines in 5 files covered. (57.14%)

527 existing lines in 19 files now uncovered.

19997 of 26966 relevant lines covered (74.16%)

67202602.73 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/math/petsc_utils/petsc_utils.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 unsigned int num_groups = uk_man_.unknowns.front().num_components;

  OpenSnPETScCall(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 (unsigned int 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);
      OpenSnPETScCall(
        MatSetValues(A_, nn, cell_idxs.data(), nn, cell_idxs.data(), cell_A.data(), ADD_VALUES));
      OpenSnPETScCall(VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES));
    } // for g
  } // for cell

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

  if (options.verbose)
  {
    MatInfo info;
    OpenSnPETScCall(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;
    OpenSnPETScCall(MatIsSymmetric(A_, 1.0e-6, &symmetry));
    if (symmetry == PETSC_FALSE)
      throw std::logic_error(fname + ":Symmetry check failed");
  }

  OpenSnPETScCall(KSPSetOperators(ksp_, A_, A_));

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

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

  OpenSnPETScCall(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 unsigned int num_groups = uk_man_.unknowns.front().num_components;

  OpenSnPETScCall(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 (unsigned int 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);
      OpenSnPETScCall(VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES));
    } // for g
  } // for cell

  OpenSnPETScCall(VecAssemblyBegin(rhs_));
  OpenSnPETScCall(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 unsigned int num_groups = uk_man_.unknowns.front().num_components;

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

  OpenSnPETScCall(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 (unsigned int 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);
      OpenSnPETScCall(VecSetValues(rhs_, nn, cell_idxs.data(), cell_rhs.data(), ADD_VALUES));
    } // for g
  } // for cell

  OpenSnPETScCall(VecRestoreArrayRead(petsc_q_vector, &q_vector));

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

Open-Sn / opensn / 22249939997

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