23371558143

Committed 20 Mar 2026 05:40PM UTC coverage: 73.918% (-0.5%) from 74.407%

Build # 23371558143

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #980 from wdhawkins/transient_fix

Transient solver should always compute t^(n+1)

Coverage Stats

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

304 existing lines in 15 files now uncovered.

20502 of 27736 relevant lines covered (73.92%)

66612898.8 hits per line

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

64.48

/framework/math/math.cc

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

#include "framework/math/math.h"
#include "framework/runtime.h"
#include "framework/logging/log.h"
#include <cassert>
#include <petscmat.h>
#include <petscksp.h>
#include <petscblaslapack.h>
#include <iostream>
#include <iomanip>
#include <limits>

namespace opensn
{

double
Factorial(const int x)
{
  double factorial_value = 1.0;
  for (int i = 2; i <= x; ++i)
    factorial_value *= i;

  return factorial_value;
}

std::pair<double, double>
OmegaToPhiThetaSafe(const Vector3& omega)
{
  // Notes: asin maps [-1,+1] to [-pi/2,+pi/2]
  //        acos maps [-1,+1] to [0,pi]
  // This mapping requires some logic for determining the azimuthal angle.
  //
  const auto omega_hat = omega.Normalized();

  double mu = omega_hat.z;
  mu = std::min(mu, 1.0);
  mu = std::max(mu, -1.0);

  double theta = acos(mu);

  // Handling omega aligned to k_hat
  if (std::fabs(omega_hat.z) < 1.0e-16)
    return {0.0, theta};

  double cos_phi = omega_hat.x / sin(theta);
  cos_phi = std::min(cos_phi, 1.0);
  cos_phi = std::max(cos_phi, -1.0);

  // Computing varphi for NE and NW quadrant
  if (omega_hat.y >= 0.0)
    return {acos(cos_phi), theta};
  // Computing varphi for SE and SW quadrant
  else
    return {2.0 * M_PI - acos(cos_phi), theta};
}

std::vector<std::vector<double>>
Transpose(const std::vector<std::vector<double>>& matrix)
{
  // Check if matrix is empty
  if (matrix.empty())
  {
    throw std::runtime_error("Cannot transpose empty matrix");
  }

  const size_t m = matrix.size();    // number of rows
  const size_t n = matrix[0].size(); // number of columns

  // Verify consistent row dimensions
  for (size_t i = 1; i < m; ++i)
  {
    if (matrix[i].size() != n)
    {
      throw std::runtime_error("Matrix must have consistent row dimensions");
    }
  }

  // Create transposed matrix with dimensions n x m
  std::vector<std::vector<double>> transposed(n, std::vector<double>(m));

  // Perform transpose: transposed[j][i] = matrix[i][j]
  for (size_t i = 0; i < m; ++i)
  {
    for (size_t j = 0; j < n; ++j)
    {
      transposed[j][i] = matrix[i][j];
    }
  }

  return transposed;
}

std::vector<std::vector<double>>
InvertMatrix(const std::vector<std::vector<double>>& matrix)
{
  const auto n = static_cast<PetscInt>(matrix.size());

  // Check if matrix is square
  if (n == 0)
  {
    throw std::runtime_error("Matrix must not be empty for inversion");
  }

  for (PetscInt i = 0; i < n; ++i)
  {
    if (std::cmp_not_equal(matrix[i].size(), n))
    {
      throw std::runtime_error("Matrix must be square for inversion");
    }
  }

  // Use SVD-based pseudo-inverse for numerical stability
  // This is more robust than LU factorization for ill-conditioned matrices

  // Create working arrays for LAPACK SVD (dgesvd)
  // A = U * S * V^T, so A^(-1) = V * S^(-1) * U^T
  std::vector<double> A_col_major(n * n);
  std::vector<double> U(n * n);
  std::vector<double> VT(n * n);
  std::vector<double> S(n);

  // Copy matrix to column-major format for LAPACK
  for (PetscInt j = 0; j < n; ++j)
  {
    for (PetscInt i = 0; i < n; ++i)
    {
      A_col_major[j * n + i] = matrix[i][j];
    }
  }

  // Call LAPACK dgesvd for SVD decomposition
  // jobu = 'A' means compute all columns of U
  // jobvt = 'A' means compute all rows of V^T
  char jobu = 'A';
  char jobvt = 'A';
  auto n_blas = static_cast<PetscBLASInt>(n);
  PetscBLASInt lwork = std::max(PetscBLASInt(1), 5 * n_blas);
  std::vector<double> work(lwork);
  PetscBLASInt info_svd = 0;

  // LAPACK SVD: dgesvd
  LAPACKgesvd_(&jobu,
               &jobvt,
               &n_blas,
               &n_blas,
               A_col_major.data(),
               &n_blas,
               S.data(),
               U.data(),
               &n_blas,
               VT.data(),
               &n_blas,
               work.data(),
               &lwork,
               &info_svd);

  if (info_svd != 0)
  {
    throw std::runtime_error("SVD decomposition failed with LAPACK error code: " +
                             std::to_string(info_svd));
  }

  // Compute condition number and determine tolerance for singular values
  double sigma_max = S[0];
  double sigma_min = S[n - 1];

  // Use machine epsilon scaled by matrix size and largest singular value as tolerance
  // This is a standard approach for pseudo-inverse computation
  const double machine_eps = std::numeric_limits<double>::epsilon();
  const double tol = static_cast<double>(n) * machine_eps * sigma_max;

  // Log condition number for debugging ill-conditioned matrices
  double condition_number = (sigma_min > tol) ? (sigma_max / sigma_min) : INFINITY;
  if (condition_number > 1.0e10)
  {
    log.Log0Warning() << "InvertMatrix: Matrix is ill-conditioned (condition number = "
                      << std::scientific << std::setprecision(3) << condition_number << ")";
  }

  // Compute pseudo-inverse: A^(-1) = V * S^(-1) * U^T
  // Since we have V^T from SVD, we need V = (V^T)^T
  // Result[i][j] = sum_k V[i][k] * (1/S[k]) * U[j][k]
  //              = sum_k VT[k][i] * (1/S[k]) * U[j][k]
  // In column-major: VT[k*n + i], U[k*n + j]

  std::vector<std::vector<double>> inverse(n, std::vector<double>(n, 0.0));

  for (PetscInt i = 0; i < n; ++i)
  {
    for (PetscInt j = 0; j < n; ++j)
    {
      double sum = 0.0;
      for (PetscInt k = 0; k < n; ++k)
      {
        // Only include contributions from singular values above tolerance
        if (S[k] > tol)
        {
          double v_ik = VT[k + i * n];
          double u_jk = U[j + k * n];
          sum += v_ik * (1.0 / S[k]) * u_jk;
        }
      }
      inverse[i][j] = sum;
    }
  }

  return inverse;
}

std::vector<std::vector<double>>
OrthogonalizeMatrixSpan(const std::vector<std::vector<double>>& matrix,
                        const std::vector<double>& weights)
{
  // Check an empty matrix
  if (matrix.empty())
  {
    throw std::runtime_error("Cannot orthogonalize empty matrix");
  }

  const size_t m = matrix.size();
  const size_t n = matrix[0].size();

  // Verify rectangular matrix
  for (size_t i = 1; i < m; ++i)
  {
    if (matrix[i].size() != n)
    {
      throw std::runtime_error("Matrix must have consistent column dimensions");
    }
  }

  // Create working copy of the matrix - we will orthogonalize its columns in place
  std::vector<std::vector<double>> Q = matrix;

  // Orthogonalize columns
  for (size_t j = 0; j < n; ++j)
  {
    // Two passes of orthogonalization for numerical stability
    for (int pass = 0; pass < 2; ++pass)
    {
      // Orthogonalize column j against all previous columns
      for (size_t i = 0; i < j; ++i)
      {
        // Compute weighted inner product <Q[:,i], Q[:,j]>
        double dot_product = 0.0;
        for (size_t row = 0; row < m; ++row)
        {
          dot_product += weights[row] * Q[row][i] * Q[row][j];
        }

        // Subtract projection: Q[:,j] -= dot_product * Q[:,i]
        for (size_t row = 0; row < m; ++row)
        {
          Q[row][j] -= dot_product * Q[row][i];
        }
      }
    }

    // Normalize column j under weighted inner product
    double norm_squared = 0.0;
    for (size_t row = 0; row < m; ++row)
    {
      norm_squared += weights[row] * Q[row][j] * Q[row][j];
    }
    double norm = std::sqrt(norm_squared);

    if (norm > 1e-14)
    {
      for (size_t row = 0; row < m; ++row)
      {
        Q[row][j] /= norm;
      }
    }
  }
  return Q;
}

void
PrintVector(const std::vector<double>& x)
{
  for (const auto& xi : x)
    std::cout << xi << ' ';
  std::cout << std::endl;
}

void
Scale(std::vector<double>& x, const double& val)
{
  for (double& xi : x)
    xi *= val;
}

void
Set(std::vector<double>& x, const double& val)
{
  for (double& xi : x)
    xi = val;
}

double
Dot(const std::vector<double>& x, const std::vector<double>& y)
{
  // Error Checks
  assert(!x.empty());
  assert(!y.empty());
  assert(x.size() == y.size());
  // Local Variables
  size_t n = x.size();
  double val = 0.0;

  for (size_t i = 0; i != n; ++i)
    val += x[i] * y[i];

  return val;
}

std::vector<double>
Mult(const std::vector<double>& x, const double& val)
{
  size_t n = x.size();
  std::vector<double> y(n);

  for (size_t i = 0; i != n; ++i)
    y[i] = val * x[i];

  return y;
}

double
L1Norm(const std::vector<double>& x)
{
  // Local Variables
  size_t n = x.size();
  double val = 0.0;

  for (size_t i = 0; i != n; ++i)
    val += std::fabs(x[i]);

  return val;
}

double
L2Norm(const std::vector<double>& x)
{
  // Local Variables
  size_t n = x.size();
  double val = 0.0;

  for (size_t i = 0; i != n; ++i)
    val += x[i] * x[i];

  return std::sqrt(val);
}

double
LInfNorm(const std::vector<double>& x)
{
  // Local Variables
  size_t n = x.size();
  double val = 0.0;

  for (size_t i = 0; i != n; ++i)
    val += std::max(std::fabs(x[i]), val);

  return val;
}

double
LpNorm(const std::vector<double>& x, const double& p)
{
  // Local Variables
  size_t n = x.size();
  double val = 0.0;

  for (size_t i = 0; i != n; ++i)
    val += std::pow(std::fabs(x[i]), p);

  return std::pow(val, 1.0 / p);
}

std::vector<double>
operator+(const std::vector<double>& a, const std::vector<double>& b)
{
  assert(a.size() == b.size());
  std::vector<double> result(a.size(), 0.0);

  for (size_t i = 0; i < a.size(); ++i)
    result[i] = a[i] + b[i];

  return result;
}

std::vector<double>
operator-(const std::vector<double>& a, const std::vector<double>& b)
{
  assert(a.size() == b.size());
  std::vector<double> result(a.size(), 0.0);

  for (size_t i = 0; i < a.size(); ++i)
    result[i] = a[i] - b[i];

  return result;
}

double
ComputeL2Change(std::vector<double>& x, std::vector<double>& y)
{
  assert(x.size() == y.size());

  double norm = 0.0;
  for (auto i = 0; i < x.size(); ++i)
  {
    double val = x[i] - y[i];
    norm += val * val;
  }

  double global_norm = 0.0;
  mpi_comm.all_reduce<double>(norm, global_norm, mpi::op::sum<double>());

  return std::sqrt(global_norm);
}

double
ComputePointwiseChange(std::vector<double>& x, std::vector<double>& y)
{
  assert(x.size() == y.size());

  double pw_change = 0.0;
  for (auto i = 0; i < x.size(); ++i)
  {
    double max = std::max(x[i], y[i]);
    double delta = std::fabs(x[i] - y[i]);
    if (max >= std::numeric_limits<double>::min())
      pw_change = std::max(delta / max, pw_change);
    else
      pw_change = std::max(delta, pw_change);
  }

  double global_pw_change = 0.0;
  mpi_comm.all_reduce<double>(pw_change, global_pw_change, mpi::op::max<double>());

  return global_pw_change;
}

} // 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/math.h"
5	#include "framework/runtime.h"
6	#include "framework/logging/log.h"
7	#include <cassert>
8	#include <petscmat.h>
9	#include <petscksp.h>
10	#include <petscblaslapack.h>
11	#include <iostream>
12	#include <iomanip>
13	#include <limits>
14
15	namespace opensn
16	{
17
18	double
19	Factorial(const int x)	30,432,784✔
20	{
21	double factorial_value = 1.0;	30,432,784✔
22	for (int i = 2; i <= x; ++i)	53,950,040✔
23	factorial_value *= i;	23,517,256✔
24
25	return factorial_value;	30,432,784✔
26	}
27
28	std::pair<double, double>
29	OmegaToPhiThetaSafe(const Vector3& omega)	839,857✔
30	{
31	// Notes: asin maps [-1,+1] to [-pi/2,+pi/2]
32	// acos maps [-1,+1] to [0,pi]
33	// This mapping requires some logic for determining the azimuthal angle.
34	//
35	const auto omega_hat = omega.Normalized();	839,857✔
36
37	double mu = omega_hat.z;	839,857✔
38	mu = std::min(mu, 1.0);	839,857✔
39	mu = std::max(mu, -1.0);	839,857✔
40
41	double theta = acos(mu);	839,857✔
42
43	// Handling omega aligned to k_hat
44	if (std::fabs(omega_hat.z) < 1.0e-16)	839,857✔
UNCOV 45	return {0.0, theta};	×
46
47	double cos_phi = omega_hat.x / sin(theta);	839,857✔
48	cos_phi = std::min(cos_phi, 1.0);	839,857✔
49	cos_phi = std::max(cos_phi, -1.0);	839,857✔
50
51	// Computing varphi for NE and NW quadrant
52	if (omega_hat.y >= 0.0)	839,857✔
53	return {acos(cos_phi), theta};	421,089✔
54	// Computing varphi for SE and SW quadrant
55	else
56	return {2.0 * M_PI - acos(cos_phi), theta};	418,768✔
57	}
58
59	std::vector<std::vector<double>>
60	Transpose(const std::vector<std::vector<double>>& matrix)	11✔
61	{
62	// Check if matrix is empty
63	if (matrix.empty())	11✔
64	{
65	throw std::runtime_error("Cannot transpose empty matrix");	1✔
66	}
67
68	const size_t m = matrix.size(); // number of rows	10✔
69	const size_t n = matrix[0].size(); // number of columns	10✔
70
71	// Verify consistent row dimensions
72	for (size_t i = 1; i < m; ++i)	626✔
73	{
74	if (matrix[i].size() != n)	616✔
75	{
UNCOV 76	throw std::runtime_error("Matrix must have consistent row dimensions");	×
77	}
78	}
79
80	// Create transposed matrix with dimensions n x m
81	std::vector<std::vector<double>> transposed(n, std::vector<double>(m));	20✔
82
83	// Perform transpose: transposed[j][i] = matrix[i][j]
84	for (size_t i = 0; i < m; ++i)	636✔
85	{
86	for (size_t j = 0; j < n; ++j)	198,152✔
87	{
88	transposed[j][i] = matrix[i][j];	197,526✔
89	}
90	}
91
92	return transposed;	10✔
93	}
94
95	std::vector<std::vector<double>>
96	InvertMatrix(const std::vector<std::vector<double>>& matrix)	8✔
97	{
98	const auto n = static_cast<PetscInt>(matrix.size());	8✔
99
100	// Check if matrix is square
101	if (n == 0)	8✔
102	{
103	throw std::runtime_error("Matrix must not be empty for inversion");	1✔
104	}
105
106	for (PetscInt i = 0; i < n; ++i)	464✔
107	{
108	if (std::cmp_not_equal(matrix[i].size(), n))	457✔
109	{
110	throw std::runtime_error("Matrix must be square for inversion");	×
111	}
112	}
113
114	// Use SVD-based pseudo-inverse for numerical stability
115	// This is more robust than LU factorization for ill-conditioned matrices
116
117	// Create working arrays for LAPACK SVD (dgesvd)
118	// A = U * S * V^T, so A^(-1) = V * S^(-1) * U^T
119	std::vector<double> A_col_major(n * n);	7✔
120	std::vector<double> U(n * n);	7✔
121	std::vector<double> VT(n * n);	7✔
122	std::vector<double> S(n);	7✔
123
124	// Copy matrix to column-major format for LAPACK
125	for (PetscInt j = 0; j < n; ++j)	464✔
126	{
127	for (PetscInt i = 0; i < n; ++i)	149,982✔
128	{
129	A_col_major[j * n + i] = matrix[i][j];	149,525✔
130	}
131	}
132
133	// Call LAPACK dgesvd for SVD decomposition
134	// jobu = 'A' means compute all columns of U
135	// jobvt = 'A' means compute all rows of V^T
136	char jobu = 'A';	7✔
137	char jobvt = 'A';	7✔
138	auto n_blas = static_cast<PetscBLASInt>(n);	7✔
139	PetscBLASInt lwork = std::max(PetscBLASInt(1), 5 * n_blas);	7✔
140	std::vector<double> work(lwork);	7✔
141	PetscBLASInt info_svd = 0;	7✔
142
143	// LAPACK SVD: dgesvd
144	LAPACKgesvd_(&jobu,	7✔
145	&jobvt,
146	&n_blas,
147	&n_blas,
148	A_col_major.data(),
149	&n_blas,
150	S.data(),
151	U.data(),
152	&n_blas,
153	VT.data(),
154	&n_blas,
155	work.data(),
156	&lwork,
157	&info_svd);
158
159	if (info_svd != 0)	7✔
160	{
UNCOV 161	throw std::runtime_error("SVD decomposition failed with LAPACK error code: " +	×
UNCOV 162	std::to_string(info_svd));	×
163	}
164
165	// Compute condition number and determine tolerance for singular values
166	double sigma_max = S[0];	7✔
167	double sigma_min = S[n - 1];	7✔
168
169	// Use machine epsilon scaled by matrix size and largest singular value as tolerance
170	// This is a standard approach for pseudo-inverse computation
171	const double machine_eps = std::numeric_limits<double>::epsilon();	7✔
172	const double tol = static_cast<double>(n) * machine_eps * sigma_max;	7✔
173
174	// Log condition number for debugging ill-conditioned matrices
175	double condition_number = (sigma_min > tol) ? (sigma_max / sigma_min) : INFINITY;	7✔
176	if (condition_number > 1.0e10)	7✔
177	{
UNCOV 178	log.Log0Warning() << "InvertMatrix: Matrix is ill-conditioned (condition number = "	×
UNCOV 179	<< std::scientific << std::setprecision(3) << condition_number << ")";	×
180	}
181
182	// Compute pseudo-inverse: A^(-1) = V * S^(-1) * U^T
183	// Since we have V^T from SVD, we need V = (V^T)^T
184	// Result[i][j] = sum_k V[i][k] * (1/S[k]) * U[j][k]
185	// = sum_k VT[k][i] * (1/S[k]) * U[j][k]
186	// In column-major: VT[kn + i], U[kn + j]
187
188	std::vector<std::vector<double>> inverse(n, std::vector<double>(n, 0.0));	14✔
189
190	for (PetscInt i = 0; i < n; ++i)	464✔
191	{
192	for (PetscInt j = 0; j < n; ++j)	149,982✔
193	{
194	double sum = 0.0;
195	for (PetscInt k = 0; k < n; ++k)	56,838,216✔
196	{
197	// Only include contributions from singular values above tolerance
198	if (S[k] > tol)	56,688,691✔
199	{
200	double v_ik = VT[k + i * n];	56,688,691✔
201	double u_jk = U[j + k * n];	56,688,691✔
202	sum += v_ik * (1.0 / S[k]) * u_jk;	56,688,691✔
203	}
204	}
205	inverse[i][j] = sum;	149,525✔
206	}
207	}
208
209	return inverse;	7✔
210	}	7✔
211
212	std::vector<std::vector<double>>
213	OrthogonalizeMatrixSpan(const std::vector<std::vector<double>>& matrix,	10✔
214	const std::vector<double>& weights)
215	{
216	// Check an empty matrix
217	if (matrix.empty())	10✔
218	{
219	throw std::runtime_error("Cannot orthogonalize empty matrix");	1✔
220	}
221
222	const size_t m = matrix.size();	9✔
223	const size_t n = matrix[0].size();	9✔
224
225	// Verify rectangular matrix
226	for (size_t i = 1; i < m; ++i)	904✔
227	{
228	if (matrix[i].size() != n)	895✔
229	{
UNCOV 230	throw std::runtime_error("Matrix must have consistent column dimensions");	×
231	}
232	}
233
234	// Create working copy of the matrix - we will orthogonalize its columns in place
235	std::vector<std::vector<double>> Q = matrix;	9✔
236
237	// Orthogonalize columns
238	for (size_t j = 0; j < n; ++j)	352✔
239	{
240	// Two passes of orthogonalization for numerical stability
241	for (int pass = 0; pass < 2; ++pass)	1,029✔
242	{
243	// Orthogonalize column j against all previous columns
244	for (size_t i = 0; i < j; ++i)	32,194✔
245	{
246	// Compute weighted inner product <Q[:,i], Q[:,j]>
247	double dot_product = 0.0;
248	for (size_t row = 0; row < m; ++row)	11,261,728✔
249	{
250	dot_product += weights[row] * Q[row][i] * Q[row][j];	11,230,220✔
251	}
252
253	// Subtract projection: Q[:,j] -= dot_product * Q[:,i]
254	for (size_t row = 0; row < m; ++row)	11,261,728✔
255	{
256	Q[row][j] -= dot_product * Q[row][i];	11,230,220✔
257	}
258	}
259	}
260
261	// Normalize column j under weighted inner product
262	double norm_squared = 0.0;
263	for (size_t row = 0; row < m; ++row)	96,357✔
264	{
265	norm_squared += weights[row] * Q[row][j] * Q[row][j];	96,014✔
266	}
267	double norm = std::sqrt(norm_squared);	343✔
268
269	if (norm > 1e-14)	343✔
270	{
271	for (size_t row = 0; row < m; ++row)	96,357✔
272	{
273	Q[row][j] /= norm;	96,014✔
274	}
275	}
276	}
277	return Q;	9✔
278	}
279
280	void
UNCOV 281	PrintVector(const std::vector<double>& x)	×
282	{
UNCOV 283	for (const auto& xi : x)	×
UNCOV 284	std::cout << xi << ' ';	×
UNCOV 285	std::cout << std::endl;	×
UNCOV 286	}	×
287
288	void
289	Scale(std::vector<double>& x, const double& val)	170✔
290	{
291	for (double& xi : x)	5,629,938✔
292	xi *= val;	5,629,768✔
293	}	170✔
294
295	void
UNCOV 296	Set(std::vector<double>& x, const double& val)	×
297	{
UNCOV 298	for (double& xi : x)	×
UNCOV 299	xi = val;	×
UNCOV 300	}	×
301
302	double
UNCOV 303	Dot(const std::vector<double>& x, const std::vector<double>& y)	×
304	{
305	// Error Checks
UNCOV 306	assert(!x.empty());	×
UNCOV 307	assert(!y.empty());	×
UNCOV 308	assert(x.size() == y.size());	×
309	// Local Variables
UNCOV 310	size_t n = x.size();	×
UNCOV 311	double val = 0.0;	×
312
UNCOV 313	for (size_t i = 0; i != n; ++i)	×
UNCOV 314	val += x[i] * y[i];	×
315
UNCOV 316	return val;	×
317	}
318
319	std::vector<double>
UNCOV 320	Mult(const std::vector<double>& x, const double& val)	×
321	{
UNCOV 322	size_t n = x.size();	×
UNCOV 323	std::vector<double> y(n);	×
324
UNCOV 325	for (size_t i = 0; i != n; ++i)	×
UNCOV 326	y[i] = val * x[i];	×
327
UNCOV 328	return y;	×
329	}
330
331	double
UNCOV 332	L1Norm(const std::vector<double>& x)	×
333	{
334	// Local Variables
UNCOV 335	size_t n = x.size();	×
UNCOV 336	double val = 0.0;	×
337
UNCOV 338	for (size_t i = 0; i != n; ++i)	×
UNCOV 339	val += std::fabs(x[i]);	×
340
UNCOV 341	return val;	×
342	}
343
344	double
UNCOV 345	L2Norm(const std::vector<double>& x)	×
346	{
347	// Local Variables
UNCOV 348	size_t n = x.size();	×
UNCOV 349	double val = 0.0;	×
350
UNCOV 351	for (size_t i = 0; i != n; ++i)	×
UNCOV 352	val += x[i] * x[i];	×
353
UNCOV 354	return std::sqrt(val);	×
355	}
356
357	double
UNCOV 358	LInfNorm(const std::vector<double>& x)	×
359	{
360	// Local Variables
UNCOV 361	size_t n = x.size();	×
UNCOV 362	double val = 0.0;	×
363
UNCOV 364	for (size_t i = 0; i != n; ++i)	×
UNCOV 365	val += std::max(std::fabs(x[i]), val);	×
366
UNCOV 367	return val;	×
368	}
369
370	double
UNCOV 371	LpNorm(const std::vector<double>& x, const double& p)	×
372	{
373	// Local Variables
UNCOV 374	size_t n = x.size();	×
UNCOV 375	double val = 0.0;	×
376
UNCOV 377	for (size_t i = 0; i != n; ++i)	×
UNCOV 378	val += std::pow(std::fabs(x[i]), p);	×
379
UNCOV 380	return std::pow(val, 1.0 / p);	×
381	}
382
383	std::vector<double>
384	operator+(const std::vector<double>& a, const std::vector<double>& b)	10,128✔
385	{
386	assert(a.size() == b.size());	10,128✔
387	std::vector<double> result(a.size(), 0.0);	10,128✔
388
389	for (size_t i = 0; i < a.size(); ++i)	36,976,528✔
390	result[i] = a[i] + b[i];	36,966,400✔
391
392	return result;	10,128✔
393	}
394
395	std::vector<double>
396	operator-(const std::vector<double>& a, const std::vector<double>& b)	2,600✔
397	{
398	assert(a.size() == b.size());	2,600✔
399	std::vector<double> result(a.size(), 0.0);	2,600✔
400
401	for (size_t i = 0; i < a.size(); ++i)	8,322,600✔
402	result[i] = a[i] - b[i];	8,320,000✔
403
404	return result;	2,600✔
405	}
406
407	double
UNCOV 408	ComputeL2Change(std::vector<double>& x, std::vector<double>& y)	×
409	{
UNCOV 410	assert(x.size() == y.size());	×
411
UNCOV 412	double norm = 0.0;	×
UNCOV 413	for (auto i = 0; i < x.size(); ++i)	×
414	{
UNCOV 415	double val = x[i] - y[i];	×
UNCOV 416	norm += val * val;	×
417	}
418
UNCOV 419	double global_norm = 0.0;	×
UNCOV 420	mpi_comm.all_reduce<double>(norm, global_norm, mpi::op::sum<double>());	×
421
UNCOV 422	return std::sqrt(global_norm);	×
423	}
424
425	double
426	ComputePointwiseChange(std::vector<double>& x, std::vector<double>& y)	148,452✔
427	{
428	assert(x.size() == y.size());	148,452✔
429
430	double pw_change = 0.0;	148,452✔
431	for (auto i = 0; i < x.size(); ++i)	102,010,868✔
432	{
433	double max = std::max(x[i], y[i]);	101,862,416✔
434	double delta = std::fabs(x[i] - y[i]);	101,862,416✔
435	if (max >= std::numeric_limits<double>::min())	101,862,416✔
436	pw_change = std::max(delta / max, pw_change);	79,080,716✔
437	else
438	pw_change = std::max(delta, pw_change);	22,781,700✔
439	}
440
441	double global_pw_change = 0.0;	148,452✔
442	mpi_comm.all_reduce<double>(pw_change, global_pw_change, mpi::op::max<double>());	148,452✔
443
444	return global_pw_change;	148,452✔
445	}
446
447	} // namespace opensn

Open-Sn / opensn / 23371558143

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