16821314466

Committed 07 Aug 2025 04:57PM UTC coverage: 74.527% (+1.0%) from 73.488%

Build # 16821314466

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #708 from wdhawkins/curvilinear_warning

Adding experimental warning to curvilinear solver

Coverage Stats

3 of 3 new or added lines in 1 file covered. (100.0%)

250 existing lines in 17 files now uncovered.

17543 of 23539 relevant lines covered (74.53%)

44793362.81 hits per line

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

79.85

/framework/math/spatial_discretization/spatial_discretization.cc

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

#include "framework/math/spatial_discretization/spatial_discretization.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/math/petsc_utils/petsc_utils.h"
#include "framework/logging/log.h"

namespace opensn
{

SpatialDiscretization::SpatialDiscretization(const std::shared_ptr<MeshContinuum> grid,
                                             SpatialDiscretizationType sdm_type)
  : UNITARY_UNKNOWN_MANAGER({std::make_pair(UnknownType::SCALAR, 0)}), grid_(grid), type_(sdm_type)
{
}

const CellMapping&
SpatialDiscretization::GetCellMapping(const Cell& cell) const
{
  constexpr std::string_view fname = "spatial_discretization::"
                                     "GetCellMapping";
  try
  {
    if (GetGrid()->IsCellLocal(cell.global_id))
      return *cell_mappings_.at(cell.local_id);
    else
      return *nb_cell_mappings_.at(cell.global_id);
  }
  catch (const std::out_of_range& oor)
  {
    throw std::out_of_range(std::string(fname) + ": Failed to obtain cell mapping.");
  }
}

SpatialDiscretizationType
SpatialDiscretization::GetType() const
{
  return type_;
}

const std::shared_ptr<MeshContinuum>
SpatialDiscretization::GetGrid() const
{
  return grid_;
}

size_t
SpatialDiscretization::GetNumLocalNodes() const
{
  return local_base_block_size_;
}

size_t
SpatialDiscretization::GetNumLocalDOFs(const UnknownManager& unknown_manager) const
{
  unsigned int N = unknown_manager.GetTotalUnknownStructureSize();

  return local_base_block_size_ * N;
}

size_t
SpatialDiscretization::GetNumGlobalNodes() const
{
  return global_base_block_size_;
}

size_t
SpatialDiscretization::GetNumGlobalDOFs(const UnknownManager& unknown_manager) const
{
  unsigned int N = unknown_manager.GetTotalUnknownStructureSize();

  return global_base_block_size_ * N;
}

size_t
SpatialDiscretization::GetNumLocalAndGhostDOFs(const UnknownManager& unknown_manager) const
{
  return GetNumLocalDOFs(unknown_manager) + GetNumGhostDOFs(unknown_manager);
}

size_t
SpatialDiscretization::GetCellNumNodes(const Cell& cell) const
{
  return GetCellMapping(cell).GetNumNodes();
}

const std::vector<Vector3>&
SpatialDiscretization::GetCellNodeLocations(const Cell& cell) const
{
  return GetCellMapping(cell).GetNodeLocations();
}

std::pair<std::set<uint32_t>, std::set<uint32_t>>
SpatialDiscretization::MakeCellInternalAndBndryNodeIDs(const Cell& cell) const
{
  const auto& cell_mapping = GetCellMapping(cell);
  const size_t num_faces = cell.faces.size();
  const size_t num_nodes = cell_mapping.GetNumNodes();

  // Determine which nodes are on the boundary
  std::set<uint32_t> boundary_nodes;
  for (size_t f = 0; f < num_faces; ++f)
  {
    if (not cell.faces[f].has_neighbor)
    {
      const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);
      for (size_t fi = 0; fi < num_face_nodes; ++fi)
        boundary_nodes.insert(cell_mapping.MapFaceNode(f, fi));
    }
  } // for f

  // Determine non-boundary nodes
  std::set<uint32_t> internal_nodes;
  for (size_t i = 0; i < num_nodes; ++i)
    if (boundary_nodes.find(i) == boundary_nodes.end())
      internal_nodes.insert(i);

  return {internal_nodes, boundary_nodes};
}

std::vector<std::vector<std::vector<int>>>
SpatialDiscretization::MakeInternalFaceNodeMappings(const double tolerance) const
{
  std::vector<std::vector<std::vector<int>>> cell_adj_mapping;
  for (const auto& cell : grid_->local_cells)
  {
    const auto& cell_mapping = this->GetCellMapping(cell);
    const auto& node_locations = cell_mapping.GetNodeLocations();
    const size_t num_faces = cell.faces.size();

    std::vector<std::vector<int>> per_face_adj_mapping;

    for (size_t f = 0; f < num_faces; ++f)
    {
      const auto& face = cell.faces[f];
      const auto num_face_nodes = cell_mapping.GetNumFaceNodes(f);
      std::vector<int> face_adj_mapping(num_face_nodes, -1);
      if (face.has_neighbor)
      {
        const auto& adj_cell = grid_->cells[face.neighbor_id];
        const auto& adj_cell_mapping = this->GetCellMapping(adj_cell);
        const auto& adj_node_locations = adj_cell_mapping.GetNodeLocations();
        const size_t adj_num_nodes = adj_cell_mapping.GetNumNodes();

        for (size_t fi = 0; fi < num_face_nodes; ++fi)
        {
          const int i = cell_mapping.MapFaceNode(f, fi);
          const auto& ivec3 = node_locations[i];

          for (size_t ai = 0; ai < adj_num_nodes; ++ai)
          {
            const auto& aivec3 = adj_node_locations[ai];
            if ((ivec3 - aivec3).NormSquare() < tolerance)
            {
              face_adj_mapping[fi] = static_cast<int>(ai);
              break;
            }
          } // for ai
          if (face_adj_mapping[fi] < 0)
            throw std::logic_error("Face node mapping failed");
        } // for fi
      } // if internal face

      per_face_adj_mapping.push_back(std::move(face_adj_mapping));
    } // for face

    cell_adj_mapping.push_back(std::move(per_face_adj_mapping));
  } // for cell

  return cell_adj_mapping;
}

void
SpatialDiscretization::CopyVectorWithUnknownScope(const std::vector<double>& from_vector,
                                                  std::vector<double>& to_vector,
                                                  const UnknownManager& from_vec_uk_structure,
                                                  const unsigned int from_vec_uk_id,
                                                  const UnknownManager& to_vec_uk_structure,
                                                  const unsigned int to_vec_uk_id) const
{
  const std::string fname = "spatial_discretization::"
                            "CopyVectorWithUnknownScope";
  const auto& ukmanF = from_vec_uk_structure;
  const auto& ukmanT = to_vec_uk_structure;
  const auto& ukidF = from_vec_uk_id;
  const auto& ukidT = to_vec_uk_id;
  try
  {
    const auto& ukA = from_vec_uk_structure.unknowns.at(from_vec_uk_id);
    const auto& ukB = to_vec_uk_structure.unknowns.at(to_vec_uk_id);

    if (ukA.num_components != ukB.num_components)
      throw std::logic_error(fname + " Unknowns do not have the "
                                     "same number of components");

    const size_t num_comps = ukA.num_components;

    for (const auto& cell : grid_->local_cells)
    {
      const auto& cell_mapping = this->GetCellMapping(cell);
      const size_t num_nodes = cell_mapping.GetNumNodes();

      for (size_t i = 0; i < num_nodes; ++i)
      {
        for (size_t c = 0; c < num_comps; ++c)
        {
          const int64_t fmap = MapDOFLocal(cell, i, ukmanF, ukidF, c);
          const int64_t imap = MapDOFLocal(cell, i, ukmanT, ukidT, c);

          to_vector[imap] = from_vector[fmap];
        } // for component c
      } // for node i
    } // for cell
  }
  catch (const std::out_of_range& oor)
  {
    throw std::out_of_range(fname + ": either from_vec_uk_id or to_vec_uk_id is "
                                    "out of range for its respective "
                                    "unknown manager.");
  }
}

void
SpatialDiscretization::LocalizePETScVector(Vec petsc_vector,
                                           std::vector<double>& local_vector,
                                           const UnknownManager& unknown_manager) const
{
  size_t num_local_dofs = GetNumLocalDOFs(unknown_manager);

  CopyVecToSTLvector(petsc_vector, local_vector, num_local_dofs);
}

void
SpatialDiscretization::LocalizePETScVectorWithGhosts(Vec petsc_vector,
                                                     std::vector<double>& local_vector,
                                                     const UnknownManager& unknown_manager) const
{
  size_t num_local_dofs = GetNumLocalAndGhostDOFs(unknown_manager);

  CopyVecToSTLvectorWithGhosts(petsc_vector, local_vector, num_local_dofs);
}

double
SpatialDiscretization::CartesianSpatialWeightFunction(const Vector3& point)
{
  return 1.0;
}

double
SpatialDiscretization::CylindricalRZSpatialWeightFunction(const Vector3& point)
{
  return 2.0 * M_PI * point[0];
}

double
SpatialDiscretization::Spherical1DSpatialWeightFunction(const Vector3& point)
{
  const double r = point[2];
  return 4.0 * M_PI * r * r;
}

SpatialDiscretization::SpatialWeightFunction
SpatialDiscretization::GetSpatialWeightingFunction() const
{
  switch (grid_->GetCoordinateSystem())
  {
    case CoordinateSystemType::CARTESIAN:
      return CartesianSpatialWeightFunction;
    case CoordinateSystemType::CYLINDRICAL:
      return CylindricalRZSpatialWeightFunction;
    case CoordinateSystemType::SPHERICAL:
      return Spherical1DSpatialWeightFunction;
    case CoordinateSystemType::UNDEFINED:
    default:
      OpenSnLogicalError("Coordinate system undefined.");
  }
}

} // namespace opensn

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "framework/math/spatial_discretization/spatial_discretization.h"
5	#include "framework/mesh/mesh_continuum/mesh_continuum.h"
6	#include "framework/math/petsc_utils/petsc_utils.h"
7	#include "framework/logging/log.h"
8
9	namespace opensn
10	{
11
12	SpatialDiscretization::SpatialDiscretization(const std::shared_ptr<MeshContinuum> grid,	339✔
13	SpatialDiscretizationType sdm_type)	339✔
14	: UNITARY_UNKNOWN_MANAGER({std::make_pair(UnknownType::SCALAR, 0)}), grid_(grid), type_(sdm_type)	339✔
15	{
16	}	339✔
17
18	const CellMapping&
19	SpatialDiscretization::GetCellMapping(const Cell& cell) const	437,452,109✔
20	{
21	constexpr std::string_view fname = "spatial_discretization::"	437,452,109✔
22	"GetCellMapping";
23	try	437,452,109✔
24	{
25	if (GetGrid()->IsCellLocal(cell.global_id))	874,904,218✔
26	return *cell_mappings_.at(cell.local_id);	436,774,867✔
27	else
28	return *nb_cell_mappings_.at(cell.global_id);	677,242✔
29	}
30	catch (const std::out_of_range& oor)	×
31	{
32	throw std::out_of_range(std::string(fname) + ": Failed to obtain cell mapping.");	×
33	}	×
34	}
35
36	SpatialDiscretizationType
37	SpatialDiscretization::GetType() const	85✔
38	{
39	return type_;	85✔
40	}
41
42	const std::shared_ptr<MeshContinuum>
43	SpatialDiscretization::GetGrid() const	437,541,984✔
44	{
45	return grid_;	437,541,984✔
46	}
47
48	size_t
49	SpatialDiscretization::GetNumLocalNodes() const	300✔
50	{
51	return local_base_block_size_;	300✔
52	}
53
54	size_t
55	SpatialDiscretization::GetNumLocalDOFs(const UnknownManager& unknown_manager) const	107,692✔
56	{
57	unsigned int N = unknown_manager.GetTotalUnknownStructureSize();	107,692✔
58
59	return local_base_block_size_ * N;	107,692✔
60	}
61
62	size_t
63	SpatialDiscretization::GetNumGlobalNodes() const	260✔
64	{
65	return global_base_block_size_;	260✔
66	}
67
68	size_t
69	SpatialDiscretization::GetNumGlobalDOFs(const UnknownManager& unknown_manager) const	72,415✔
70	{
71	unsigned int N = unknown_manager.GetTotalUnknownStructureSize();	72,415✔
72
73	return global_base_block_size_ * N;	72,415✔
74	}
75
76	size_t
77	SpatialDiscretization::GetNumLocalAndGhostDOFs(const UnknownManager& unknown_manager) const	133✔
78	{
79	return GetNumLocalDOFs(unknown_manager) + GetNumGhostDOFs(unknown_manager);	133✔
80	}
81
82	size_t
83	SpatialDiscretization::GetCellNumNodes(const Cell& cell) const	1,768,448✔
84	{
85	return GetCellMapping(cell).GetNumNodes();	1,768,448✔
86	}
87
88	const std::vector<Vector3>&
89	SpatialDiscretization::GetCellNodeLocations(const Cell& cell) const	1,099,126✔
90	{
91	return GetCellMapping(cell).GetNodeLocations();	1,099,126✔
92	}
93
94	std::pair<std::set<uint32_t>, std::set<uint32_t>>
95	SpatialDiscretization::MakeCellInternalAndBndryNodeIDs(const Cell& cell) const	141,912✔
96	{
97	const auto& cell_mapping = GetCellMapping(cell);	141,912✔
98	const size_t num_faces = cell.faces.size();	141,912✔
99	const size_t num_nodes = cell_mapping.GetNumNodes();	141,912✔
100
101	// Determine which nodes are on the boundary
102	std::set<uint32_t> boundary_nodes;	141,912✔
103	for (size_t f = 0; f < num_faces; ++f)	892,896✔
104	{
105	if (not cell.faces[f].has_neighbor)	750,984✔
106	{
107	const size_t num_face_nodes = cell_mapping.GetNumFaceNodes(f);	21,540✔
108	for (size_t fi = 0; fi < num_face_nodes; ++fi)	98,044✔
109	boundary_nodes.insert(cell_mapping.MapFaceNode(f, fi));	76,504✔
110	}
111	} // for f
112
113	// Determine non-boundary nodes
114	std::set<uint32_t> internal_nodes;	141,912✔
115	for (size_t i = 0; i < num_nodes; ++i)	1,077,872✔
116	if (boundary_nodes.find(i) == boundary_nodes.end())	935,960✔
117	internal_nodes.insert(i);	862,022✔
118
119	return {internal_nodes, boundary_nodes};	141,912✔
120	}	141,912✔
121
122	std::vector<std::vector<std::vector<int>>>
123	SpatialDiscretization::MakeInternalFaceNodeMappings(const double tolerance) const	1✔
124	{
125	std::vector<std::vector<std::vector<int>>> cell_adj_mapping;	1✔
126	for (const auto& cell : grid_->local_cells)	626✔
127	{
128	const auto& cell_mapping = this->GetCellMapping(cell);	625✔
129	const auto& node_locations = cell_mapping.GetNodeLocations();	625✔
130	const size_t num_faces = cell.faces.size();	625✔
131
132	std::vector<std::vector<int>> per_face_adj_mapping;	625✔
133
134	for (size_t f = 0; f < num_faces; ++f)	3,125✔
135	{
136	const auto& face = cell.faces[f];	2,500✔
137	const auto num_face_nodes = cell_mapping.GetNumFaceNodes(f);	2,500✔
138	std::vector<int> face_adj_mapping(num_face_nodes, -1);	2,500✔
139	if (face.has_neighbor)	2,500✔
140	{
141	const auto& adj_cell = grid_->cells[face.neighbor_id];	2,400✔
142	const auto& adj_cell_mapping = this->GetCellMapping(adj_cell);	2,400✔
143	const auto& adj_node_locations = adj_cell_mapping.GetNodeLocations();	2,400✔
144	const size_t adj_num_nodes = adj_cell_mapping.GetNumNodes();	2,400✔
145
146	for (size_t fi = 0; fi < num_face_nodes; ++fi)	7,200✔
147	{
148	const int i = cell_mapping.MapFaceNode(f, fi);	4,800✔
149	const auto& ivec3 = node_locations[i];	4,800✔
150
151	for (size_t ai = 0; ai < adj_num_nodes; ++ai)	12,000✔
152	{
153	const auto& aivec3 = adj_node_locations[ai];	12,000✔
154	if ((ivec3 - aivec3).NormSquare() < tolerance)	12,000✔
155	{
156	face_adj_mapping[fi] = static_cast<int>(ai);	4,800✔
157	break;	4,800✔
158	}
159	} // for ai
160	if (face_adj_mapping[fi] < 0)	4,800✔
161	throw std::logic_error("Face node mapping failed");	×
162	} // for fi
163	} // if internal face
164
165	per_face_adj_mapping.push_back(std::move(face_adj_mapping));	2,500✔
166	} // for face	2,500✔
167
168	cell_adj_mapping.push_back(std::move(per_face_adj_mapping));	625✔
169	} // for cell	625✔
170
171	return cell_adj_mapping;	1✔
172	}	×
173
174	void
175	SpatialDiscretization::CopyVectorWithUnknownScope(const std::vector<double>& from_vector,	6✔
176	std::vector<double>& to_vector,
177	const UnknownManager& from_vec_uk_structure,
178	const unsigned int from_vec_uk_id,
179	const UnknownManager& to_vec_uk_structure,
180	const unsigned int to_vec_uk_id) const
181	{
182	const std::string fname = "spatial_discretization::"	6✔
183	"CopyVectorWithUnknownScope";	6✔
184	const auto& ukmanF = from_vec_uk_structure;	6✔
185	const auto& ukmanT = to_vec_uk_structure;	6✔
186	const auto& ukidF = from_vec_uk_id;	6✔
187	const auto& ukidT = to_vec_uk_id;	6✔
188	try	6✔
189	{
190	const auto& ukA = from_vec_uk_structure.unknowns.at(from_vec_uk_id);	6✔
191	const auto& ukB = to_vec_uk_structure.unknowns.at(to_vec_uk_id);	6✔
192
193	if (ukA.num_components != ukB.num_components)	6✔
194	throw std::logic_error(fname + " Unknowns do not have the "	×
195	"same number of components");	×
196
197	const size_t num_comps = ukA.num_components;	6✔
198
199	for (const auto& cell : grid_->local_cells)	2,727✔
200	{
201	const auto& cell_mapping = this->GetCellMapping(cell);	2,721✔
202	const size_t num_nodes = cell_mapping.GetNumNodes();	2,721✔
203
204	for (size_t i = 0; i < num_nodes; ++i)	13,305✔
205	{
206	for (size_t c = 0; c < num_comps; ++c)	213,068✔
207	{
208	const int64_t fmap = MapDOFLocal(cell, i, ukmanF, ukidF, c);	202,484✔
209	const int64_t imap = MapDOFLocal(cell, i, ukmanT, ukidT, c);	202,484✔
210
211	to_vector[imap] = from_vector[fmap];	202,484✔
212	} // for component c
213	} // for node i
214	} // for cell
215	}
216	catch (const std::out_of_range& oor)	×
217	{
218	throw std::out_of_range(fname + ": either from_vec_uk_id or to_vec_uk_id is "	×
219	"out of range for its respective "
220	"unknown manager.");	×
221	}	×
222	}	6✔
223
224	void
225	SpatialDiscretization::LocalizePETScVector(Vec petsc_vector,	9,033✔
226	std::vector<double>& local_vector,
227	const UnknownManager& unknown_manager) const
228	{
229	size_t num_local_dofs = GetNumLocalDOFs(unknown_manager);	9,033✔
230
231	CopyVecToSTLvector(petsc_vector, local_vector, num_local_dofs);	9,033✔
232	}	9,033✔
233
234	void
235	SpatialDiscretization::LocalizePETScVectorWithGhosts(Vec petsc_vector,	130✔
236	std::vector<double>& local_vector,
237	const UnknownManager& unknown_manager) const
238	{
239	size_t num_local_dofs = GetNumLocalAndGhostDOFs(unknown_manager);	130✔
240
241	CopyVecToSTLvectorWithGhosts(petsc_vector, local_vector, num_local_dofs);	130✔
242	}	130✔
243
244	double
UNCOV 245	SpatialDiscretization::CartesianSpatialWeightFunction(const Vector3& point)	×
246	{
UNCOV 247	return 1.0;	×
248	}
249
250	double
251	SpatialDiscretization::CylindricalRZSpatialWeightFunction(const Vector3& point)	3,840,000✔
252	{
253	return 2.0 * M_PI * point[0];	3,840,000✔
254	}
255
256	double
257	SpatialDiscretization::Spherical1DSpatialWeightFunction(const Vector3& point)	×
258	{
259	const double r = point[2];	×
260	return 4.0 * M_PI * r * r;	×
261	}
262
263	SpatialDiscretization::SpatialWeightFunction
UNCOV 264	SpatialDiscretization::GetSpatialWeightingFunction() const	×
265	{
UNCOV 266	switch (grid_->GetCoordinateSystem())	×
267	{
UNCOV 268	case CoordinateSystemType::CARTESIAN:	×
UNCOV 269	return CartesianSpatialWeightFunction;	×
270	case CoordinateSystemType::CYLINDRICAL:	×
271	return CylindricalRZSpatialWeightFunction;	×
272	case CoordinateSystemType::SPHERICAL:	×
273	return Spherical1DSpatialWeightFunction;	×
274	case CoordinateSystemType::UNDEFINED:	×
275	default:	×
276	OpenSnLogicalError("Coordinate system undefined.");	×
277	}
278	}
279
280	} // namespace opensn

Open-Sn / opensn / 16821314466

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