18300593117

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

Build # 18300593117

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #759 from wdhawkins/performance

Sweep performance optimizations

Coverage Stats

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

334 existing lines in 80 files now uncovered.

17788 of 23761 relevant lines covered (74.86%)

61852783.95 hits per line

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

55.56

/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 size_t 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) = 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

      MatSetValues(
        A_, num_nodes, cell_idxs.data(), num_nodes, cell_idxs.data(), cell_A.data(), ADD_VALUES);
      VecSetValues(rhs_, num_nodes, 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;
  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 size_t 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) = 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
      VecSetValues(rhs_, num_nodes, 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;
  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 size_t 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) = 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
      VecSetValues(rhs_, num_nodes, 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),	1✔
26	sdm,
27	uk_man,
28	std::move(bcs),	1✔
29	std::move(map_mat_id_2_xs),	1✔
30	unit_cell_matrices,
31	suppress_bcs,
32	true,
33	verbose)	3✔
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::"
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";	1✔
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 size_t 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) = sdm_.MapDOF(cell, i, uk_man_, 0, g);	800✔
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;	400✔
116	for (size_t j = 0; j < num_nodes; ++j)	2,000✔
117	{
118
119	const double entry_aij =
120	Dg * intV_gradshapeI_gradshapeJ(i, j) + sigr_g * intV_shapeI_shapeJ(i, j);	3,200✔
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];	3,200✔
131	} // for j
132
133	cell_rhs(i) = entry_rhs_i;	800✔
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	{
UNCOV 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);	320✔
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	MatSetValues(	100✔
203	A_, num_nodes, cell_idxs.data(), num_nodes, cell_idxs.data(), cell_A.data(), ADD_VALUES);	100✔
204	VecSetValues(rhs_, num_nodes, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);	100✔
205	} // for g	100✔
206	} // for cell	100✔
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;
216	MatGetInfo(A_, MAT_GLOBAL_SUM, &info);	1✔
217
218	log.Log() << "Number of mallocs used = " << info.mallocs	2✔
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;	1✔
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";	1✔
236
237	PC pc;
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::"
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";	1✔
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 size_t 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) = sdm_.MapDOF(cell, i, uk_man_, 0, g);	800✔
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;	400✔
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];	3,200✔
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	{
UNCOV 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	VecSetValues(rhs_, num_nodes, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);	100✔
388	} // for g	100✔
389	} // for cell	100✔
390
391	VecAssemblyBegin(rhs_);	1✔
392	VecAssemblyEnd(rhs_);	1✔
393
394	if (options.verbose)	1✔
395	log.Log() << program_timer.GetTimeString() << " Assembly completed";	1✔
396	}	1✔
397
398	void
399	DiffusionPWLCSolver::Assemble_b(Vec petsc_q_vector)	×
400	{
401	const std::string fname = "acceleration::DiffusionMIPSolver::"
402	"Assemble_b";	×
403	if (A_ == nullptr or rhs_ == nullptr or ksp_ == nullptr)	×
404	throw std::logic_error(fname + ": Some or all PETSc elements are null. "	×
405	"Check that Initialize has been called.");	×
406	if (options.verbose)	×
407	log.Log() << program_timer.GetTimeString() << " Starting assembly";	×
408
409	const size_t num_groups = uk_man_.unknowns.front().num_components;	×
410
411	const double* q_vector;
412	VecGetArrayRead(petsc_q_vector, &q_vector);	×
413
414	VecSet(rhs_, 0.0);	×
415	for (const auto& cell : grid_->local_cells)	×
416	{
417	const size_t num_faces = cell.faces.size();	×
418	const auto& cell_mapping = sdm_.GetCellMapping(cell);	×
419	const size_t num_nodes = cell_mapping.GetNumNodes();	×
420	const auto cc_nodes = cell_mapping.GetNodeLocations();	×
421	const auto& unit_cell_matrices = unit_cell_matrices_[cell.local_id];	×
422
423	const auto& intV_shapeI_shapeJ = unit_cell_matrices.intV_shapeI_shapeJ;	×
424	const auto& intV_shapeI = unit_cell_matrices.intV_shapeI;	×
425
426	// Mark dirichlet nodes
427	std::vector<bool> node_is_dirichlet(num_nodes, false);	×
428	for (size_t f = 0; f < num_faces; ++f)	×
429	{
430	const auto& face = cell.faces[f];	×
431	if (not face.has_neighbor and not suppress_bcs_)	×
432	{
433	BoundaryCondition bc;	×
434	if (bcs_.count(face.neighbor_id) > 0)	×
435	bc = bcs_.at(face.neighbor_id);	×
436
437	if (bc.type != BCType::DIRICHLET)	×
438	continue;	×
439
440	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	×
441	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
442	node_is_dirichlet[cell_mapping.MapFaceNode(f, fi)] = true;	×
443	}
444	}
445
446	Vector<double> cell_rhs(num_nodes);	×
447	Vector<PetscInt> cell_idxs(num_nodes);	×
448
449	for (size_t g = 0; g < num_groups; ++g)	×
450	{
451	cell_rhs.Set(0.);	×
452	for (size_t i = 0; i < num_nodes; ++i)	×
453	cell_idxs(i) = sdm_.MapDOF(cell, i, uk_man_, 0, g);	×
454
455	// Get coefficient and nodal src
456	std::vector<double> qg(num_nodes, 0.0);	×
457	for (size_t j = 0; j < num_nodes; ++j)	×
458	qg[j] = q_vector[sdm_.MapDOFLocal(cell, j, uk_man_, 0, g)];	×
459
460	// Assemble continuous terms
461	for (size_t i = 0; i < num_nodes; ++i)	×
462	{
463	if (node_is_dirichlet[i])	×
464	continue;	×
UNCOV 465	double entry_rhs_i = 0.0;	×
466	for (size_t j = 0; j < num_nodes; ++j)	×
467	entry_rhs_i += intV_shapeI_shapeJ(i, j) * qg[j];	×
468
469	cell_rhs(i) += entry_rhs_i;	×
470	} // for i
471
472	// Assemble face terms
473	for (size_t f = 0; f < num_faces; ++f)	×
474	{
475	const auto& face = cell.faces[f];	×
476	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	×
477
478	const auto& intS_shapeI = unit_cell_matrices.intS_shapeI[f];	×
479
480	if (not face.has_neighbor and not suppress_bcs_)	×
481	{
482	BoundaryCondition bc;	×
483	if (bcs_.count(face.neighbor_id) > 0)	×
484	bc = bcs_.at(face.neighbor_id);	×
485
486	if (bc.type == BCType::DIRICHLET)	×
487	{
UNCOV 488	const double bc_value = bc.values[0];	×
489
490	// Assembly penalty terms
491	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
492	{
493	const int i = cell_mapping.MapFaceNode(f, fi);	×
494
495	cell_rhs(i) += bc_value * intV_shapeI(i);	×
496	} // for fi
497
498	} // Dirichlet BC
499	else if (bc.type == BCType::ROBIN)	×
500	{
501	const double bval = bc.values[1];	×
502	const double fval = bc.values[2];	×
503
504	if (std::fabs(bval) < 1.0e-12)	×
505	continue; // a and f assumed zero	×
506
507	for (size_t fi = 0; fi < num_face_nodes; ++fi)	×
508	{
509	const int i = cell_mapping.MapFaceNode(f, fi);	×
510
511	if (std::fabs(fval) >= 1.0e-12)	×
512	{
513	const double rhs_val = (fval / bval) * intS_shapeI(i);	×
514
515	cell_rhs(i) += rhs_val;	×
516	} // if f nonzero
517	} // for fi
518	} // Robin BC
519	} // boundary face
520	} // for face
521	VecSetValues(rhs_, num_nodes, cell_idxs.data(), cell_rhs.data(), ADD_VALUES);	×
522	} // for g	×
523	} // for cell	×
524
525	VecRestoreArrayRead(petsc_q_vector, &q_vector);	×
526
527	VecAssemblyBegin(rhs_);	×
528	VecAssemblyEnd(rhs_);	×
529
530	if (options.verbose)	×
531	log.Log() << program_timer.GetTimeString() << " Assembly completed";	×
532	}	×
533
534	} // namespace opensn

Open-Sn / opensn / 18300593117

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

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