29167401541

Committed 11 Jul 2026 08:38PM UTC coverage: 80.762% (-0.5%) from 81.295%

Build # 29167401541

Build Type

Pull #4010

github

Committed by

web-flow

Commit Message

Merge 612540872 into e783e0147

Pull Request Pull Request #4010: Implementation of surfaces with custom implicit function

Coverage Stats

18489 of 27369 branches covered (67.55%)

Branch coverage included in aggregate %.

1041 of 1330 new or added lines in 13 files covered. (78.27%)

76 existing lines in 3 files now uncovered.

60486 of 70418 relevant lines covered (85.9%)

53343422.58 hits per line

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

45.75

/src/implicit_solvers.cpp

#include "openmc/implicit_solvers.h"

#include <fmt/core.h>

#include "openmc/constants.h"
#include "openmc/error.h"

namespace openmc {

std::unique_ptr<ImplicitSolver> ImplicitSolver::create(
  const std::string& name, int max_iter, double atol, double ftol)
{
  if (name == "naive")
    return std::make_unique<NaiveLipschitz>(atol, ftol, max_iter);
  if (name == "fast")
    return std::make_unique<FastLipschitz>(atol, ftol, max_iter);

  fatal_error("make_implicit_solver: unknown solver name '" + name +
              "'. "
              "Valid options are: 'naive', 'simple', 'fast'.");
  return nullptr; // unreachable — suppresses compiler warning
}

double NaiveLipschitz::solve(const Implicit& function, Position r, Direction u,
  double t0, double t1, double isovalue) const
{
  const double L = function.compute_lipschitz(r, u, t0, t1);

  // Constant function — root only if already within tolerance at t0.
  if (L == 0.0) {
    double f0 = function.along_ray(t0, r, u) - isovalue;
    return (std::abs(f0) < ftol_) ? t0 : INFTY;
  }

  double t_curr = t0;
  double f_curr = function.along_ray(t_curr, r, u) - isovalue;

  for (int i = 0; i < max_iter_; ++i) {
    // Safe Lipschitz step: f can only change by L*step over this distance,
    // so no root is possible before t_curr + |f_curr| / L.
    const double step = std::abs(f_curr) / L;
    const double t_next = t_curr + step;

    // Stepped past the end of the interval — check endpoint.
    if (t_next >= t1) {
      double f1 = function.along_ray(t1, r, u) - isovalue;
      if (f_curr * f1 < 0.0)
        return bisect(function, r, u, t_curr, t1, f_curr, f1, isovalue);
      return INFTY;
    }

    const double f_next = function.along_ray(t_next, r, u) - isovalue;

    if (std::abs(f_next) < ftol_)
      return t_next;

    // Sign change detected — root is bracketed in [t_curr, t_next].
    if (f_curr * f_next < 0.0)
      return bisect(function, r, u, t_curr, t_next, f_curr, f_next, isovalue);

    t_curr = t_next;
    f_curr = f_next;
  }

  // Max iterations reached — skip this surface.
  warning(fmt::format("NaiveLipschitz reached max iterations ({})."
                      " The surface will be skipped. Results could be biased.",
    max_iter_));
  return INFTY;
}

double NaiveLipschitz::bisect(const Implicit& function, Position r, Direction u,
  double ta, double tb, double fa, double fb, double isovalue) const
{
  // Invariant: fa and fb have opposite signs, root is in (ta, tb).
  // tb is always on the destination side of the crossing.
  // We return tb rather than the midpoint so that the returned position
  // is unambiguously on the destination side — this prevents sense()
  // from placing the particle back in the cell it just left.
  while ((tb - ta) > atol_) {
    const double tm = 0.5 * (ta + tb);
    const double fm = function.along_ray(tm, r, u) - isovalue;
    if (fa * fm < 0.0) {
      tb = tm;
      fb = fm;
    } else {
      ta = tm;
      fa = fm;
    }
  }
  return tb;
}

double FastLipschitz::solve(const Implicit& function, Position r, Direction u,
  double t0, double t1, double isovalue) const
{
  // Stack of intervals to explore.  Stored as (t0, t1) pairs.
  // Because the stack is LIFO, pushing the right half first ensures
  // the left half (smaller t) is always explored first, so the
  // solver finds the smallest root.
  std::vector<std::pair<double, double>> stack;
  stack.reserve(64);
  stack.push_back({t0, t1});

  int n_iter = 0;

  while (!stack.empty()) {
    auto [t0, t1] = stack.back();
    stack.pop_back();

    // Interval has collapsed to within tolerance — root is at t0.
    if (t1 - t0 < atol_)
      return t1;

    // Compute tight Lipschitz bound and function range over [t0, t1].
    const double L = function.compute_lipschitz(r, u, t0, t1);

    if (L == 0.0) {
      // Constant on this interval — root only if already on surface.
      double f0 = function.along_ray(t0, r, u) - isovalue;
      if (std::abs(f0) < atol_)
        return t0;
      continue;
    }

    auto [fmin_raw, fmax_raw] = function.compute_f_min_max(r, u, t0, t1);
    const double fmin = fmin_raw - isovalue;
    const double fmax = fmax_raw - isovalue;

    // Both bounds have the same sign — no root in this interval.
    if (fmin * fmax > 0.0)
      continue;

    // Endpoint values (needed to compute safe sub-interval).
    const double f0 = function.along_ray(t0, r, u) - isovalue;
    const double f1 = function.along_ray(t1, r, u) - isovalue;

    // Safe start: root cannot exist before tSafeStart.
    // Derivation: |f(t)| >= |f0| - L*(t-t0), so f cannot cross zero
    // until |f0| - L*(t-t0) <= 0  =>  t >= t0 + |f0|/L.
    // Subtract L*atol_ to account for the tolerance band.
    const double tSafeStart = t0 + std::max(0.0, std::abs(f0) - L * atol_) / L;
    if (tSafeStart - t0 < atol_)
      return tSafeStart; // root is within atol of t0

    // Safe end: root cannot exist after tSafeEnd (symmetric argument).
    const double tSafeEnd = t1 - std::max(0.0, std::abs(f1) - L * atol_) / L;
    if (t1 - tSafeEnd < atol_)
      return t1; // root is within atol of t1

    // Safe interval has collapsed — root is at tSafeStart.
    if (tSafeEnd - tSafeStart < atol_)
      return tSafeEnd;

    if (++n_iter >= max_iter_) {
      warning(
        fmt::format("FastLipschitz reached max iterations ({})."
                    " The surface will be skipped. Results could be biased.",
          max_iter_));
      return INFTY;
    }

    // Split the safe interval and push both halves.
    // Right is pushed first so that left is popped and explored first.
    const double tSplit = 0.5 * (tSafeStart + tSafeEnd);
    stack.push_back({tSplit, tSafeEnd});   // right — explored second
    stack.push_back({tSafeStart, tSplit}); // left  — explored first
  }

  return INFTY;
}

} // namespace openmc

1	#include "openmc/implicit_solvers.h"
2
3	#include <fmt/core.h>
4
5	#include "openmc/constants.h"
6	#include "openmc/error.h"
7
8	namespace openmc {
9
10	std::unique_ptr<ImplicitSolver> ImplicitSolver::create(	120✔
11	const std::string& name, int max_iter, double atol, double ftol)
12	{
13	if (name == "naive")	120!
NEW 14	return std::make_unique<NaiveLipschitz>(atol, ftol, max_iter);	×
15	if (name == "fast")	120!
16	return std::make_unique<FastLipschitz>(atol, ftol, max_iter);	120✔
17
NEW 18	fatal_error("make_implicit_solver: unknown solver name '" + name +	×
19	"'. "
20	"Valid options are: 'naive', 'simple', 'fast'.");
21	return nullptr; // unreachable — suppresses compiler warning
22	}
23
NEW 24	double NaiveLipschitz::solve(const Implicit& function, Position r, Direction u,	×
25	double t0, double t1, double isovalue) const
26	{
NEW 27	const double L = function.compute_lipschitz(r, u, t0, t1);	×
28
29	// Constant function — root only if already within tolerance at t0.
NEW 30	if (L == 0.0) {	×
NEW 31	double f0 = function.along_ray(t0, r, u) - isovalue;	×
NEW 32	return (std::abs(f0) < ftol_) ? t0 : INFTY;	×
33	}
34
NEW 35	double t_curr = t0;	×
NEW 36	double f_curr = function.along_ray(t_curr, r, u) - isovalue;	×
37
NEW 38	for (int i = 0; i < max_iter_; ++i) {	×
39	// Safe Lipschitz step: f can only change by L*step over this distance,
40	// so no root is possible before t_curr + \|f_curr\| / L.
NEW 41	const double step = std::abs(f_curr) / L;	×
NEW 42	const double t_next = t_curr + step;	×
43
44	// Stepped past the end of the interval — check endpoint.
NEW 45	if (t_next >= t1) {	×
NEW 46	double f1 = function.along_ray(t1, r, u) - isovalue;	×
NEW 47	if (f_curr * f1 < 0.0)	×
NEW 48	return bisect(function, r, u, t_curr, t1, f_curr, f1, isovalue);	×
49	return INFTY;
50	}
51
NEW 52	const double f_next = function.along_ray(t_next, r, u) - isovalue;	×
53
NEW 54	if (std::abs(f_next) < ftol_)	×
NEW 55	return t_next;	×
56
57	// Sign change detected — root is bracketed in [t_curr, t_next].
NEW 58	if (f_curr * f_next < 0.0)	×
NEW 59	return bisect(function, r, u, t_curr, t_next, f_curr, f_next, isovalue);	×
60
NEW 61	t_curr = t_next;	×
NEW 62	f_curr = f_next;	×
63	}
64
65	// Max iterations reached — skip this surface.
NEW 66	warning(fmt::format("NaiveLipschitz reached max iterations ({})."	×
67	" The surface will be skipped. Results could be biased.",
NEW 68	max_iter_));	×
NEW 69	return INFTY;	×
70	}
71
NEW 72	double NaiveLipschitz::bisect(const Implicit& function, Position r, Direction u,	×
73	double ta, double tb, double fa, double fb, double isovalue) const
74	{
75	// Invariant: fa and fb have opposite signs, root is in (ta, tb).
76	// tb is always on the destination side of the crossing.
77	// We return tb rather than the midpoint so that the returned position
78	// is unambiguously on the destination side — this prevents sense()
79	// from placing the particle back in the cell it just left.
NEW 80	while ((tb - ta) > atol_) {	×
NEW 81	const double tm = 0.5 * (ta + tb);	×
NEW 82	const double fm = function.along_ray(tm, r, u) - isovalue;	×
NEW 83	if (fa * fm < 0.0) {	×
84	tb = tm;
NEW 85	fb = fm;	×
86	} else {
NEW 87	ta = tm;	×
NEW 88	fa = fm;	×
89	}
90	}
NEW 91	return tb;	×
92	}
93
94	double FastLipschitz::solve(const Implicit& function, Position r, Direction u,	29,493,563✔
95	double t0, double t1, double isovalue) const
96	{
97	// Stack of intervals to explore. Stored as (t0, t1) pairs.
98	// Because the stack is LIFO, pushing the right half first ensures
99	// the left half (smaller t) is always explored first, so the
100	// solver finds the smallest root.
101	std::vector<std::pair<double, double>> stack;	29,493,563✔
102	stack.reserve(64);	29,493,563✔
103	stack.push_back({t0, t1});	29,493,563✔
104
105	int n_iter = 0;	29,493,563✔
106
107	while (!stack.empty()) {	1,257,803,360✔
108	auto [t0, t1] = stack.back();	1,255,925,132✔
109	stack.pop_back();	1,255,925,132✔
110
111	// Interval has collapsed to within tolerance — root is at t0.
112	if (t1 - t0 < atol_)	1,255,925,132✔
113	return t1;	27,615,335✔
114
115	// Compute tight Lipschitz bound and function range over [t0, t1].
116	const double L = function.compute_lipschitz(r, u, t0, t1);	1,255,920,732✔
117
118	if (L == 0.0) {	1,255,920,732!
119	// Constant on this interval — root only if already on surface.
NEW 120	double f0 = function.along_ray(t0, r, u) - isovalue;	×
NEW 121	if (std::abs(f0) < atol_)	×
NEW 122	return t0;	×
123	continue;	478,267,537✔
NEW 124	}	×
125
126	auto [fmin_raw, fmax_raw] = function.compute_f_min_max(r, u, t0, t1);	1,255,920,732✔
127	const double fmin = fmin_raw - isovalue;	1,255,920,732✔
128	const double fmax = fmax_raw - isovalue;	1,255,920,732✔
129
130	// Both bounds have the same sign — no root in this interval.
131	if (fmin * fmax > 0.0)	1,255,920,732✔
132	continue;	478,267,537✔
133
134	// Endpoint values (needed to compute safe sub-interval).
135	const double f0 = function.along_ray(t0, r, u) - isovalue;	777,653,195✔
136	const double f1 = function.along_ray(t1, r, u) - isovalue;	777,653,195✔
137
138	// Safe start: root cannot exist before tSafeStart.
139	// Derivation: \|f(t)\| >= \|f0\| - L*(t-t0), so f cannot cross zero
140	// until \|f0\| - L*(t-t0) <= 0 => t >= t0 + \|f0\|/L.
141	// Subtract L*atol_ to account for the tolerance band.
142	const double tSafeStart = t0 + std::max(0.0, std::abs(f0) - L * atol_) / L;	777,653,195✔
143	if (tSafeStart - t0 < atol_)	777,653,195✔
144	return tSafeStart; // root is within atol of t0	11,357,313✔
145
146	// Safe end: root cannot exist after tSafeEnd (symmetric argument).
147	const double tSafeEnd = t1 - std::max(0.0, std::abs(f1) - L * atol_) / L;	766,295,882✔
148	if (t1 - tSafeEnd < atol_)	766,295,882✔
149	return t1; // root is within atol of t1	16,253,611✔
150
151	// Safe interval has collapsed — root is at tSafeStart.
152	if (tSafeEnd - tSafeStart < atol_)	750,042,271✔
153	return tSafeEnd;	11✔
154
155	if (++n_iter >= max_iter_) {	750,042,260!
NEW 156	warning(	×
NEW 157	fmt::format("FastLipschitz reached max iterations ({})."	×
158	" The surface will be skipped. Results could be biased.",
NEW 159	max_iter_));	×
NEW 160	return INFTY;	×
161	}
162
163	// Split the safe interval and push both halves.
164	// Right is pushed first so that left is popped and explored first.
165	const double tSplit = 0.5 * (tSafeStart + tSafeEnd);	750,042,260✔
166	stack.push_back({tSplit, tSafeEnd}); // right — explored second	750,042,260✔
167	stack.push_back({tSafeStart, tSplit}); // left — explored first	750,042,260✔
168	}
169
170	return INFTY;
171	}	29,493,563✔
172
173	} // namespace openmc

openmc-dev / openmc / 29167401541

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