Skip to content

Reject a profile that its evaluator cannot evaluate #816

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CharlesCNorton wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Reject a profile that its evaluator cannot evaluate #816

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1b8ec3d
Reject a spline profile with too few knots to evaluate
CharlesCNorton Sep 7, 2026
c152e8b
Reject a closed-form profile short of coefficients and a spline with …
CharlesCNorton Sep 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 53 additions & 14 deletions src/vmecpp/cpp/vmecpp/common/vmec_indata/vmec_indata.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,17 +44,19 @@ std::string ProfileTypeName(vmecpp::ProfileType profile_type) {
}

// Checks that `type_name` names a profile parameterization that may be used for
// `profile_type`, and that a spline parameterization was given its knots. An
// unrecognized name otherwise reaches the solver as a zero profile, which
// converges to a silently wrong equilibrium.
// `profile_type`, that a spline parameterization was given its knots, and that
// a closed-form parameterization was given the coefficients it needs. An
// unrecognized name or a profile short of its data otherwise reaches the
// solver as a zero profile, which converges to a silently wrong equilibrium.
//
// The polynomial coefficient arrays are deliberately not required to be
// non-empty: they are zero-padded on read, so an empty array is a valid way to
// specify a zero profile.
// The polynomial coefficient arrays are zero-padded on read, so an empty
// array is a valid way to specify a zero power series; the closed forms that
// raise to a coefficient or divide by one need their full set.
absl::Status CheckProfile(const std::string& type_key,
const std::string& type_name,
vmecpp::ProfileType profile_type,
const std::string& aux_key,
const std::string& coefficient_key,
const Eigen::VectorXd& coefficients,
const Eigen::VectorXd& aux_s,
const Eigen::VectorXd& aux_f) {
const vmecpp::ProfileParameterizationData* const parameterization =
Expand All @@ -73,18 +75,55 @@ absl::Status CheckProfile(const std::string& type_key,
type_key, type_name, ProfileTypeName(profile_type)));
}

const int minimum_coefficients = parameterization->MinimumCoefficients();
if (coefficients.size() < minimum_coefficients) {
return absl::InvalidArgumentError(absl::StrFormat(
"'%s' is '%s', which needs at least %d coefficients, but '%s' has "
"%d\n",
type_key, type_name, minimum_coefficients, coefficient_key,
coefficients.size()));
}

if (parameterization->NeedsSplineData()) {
if (aux_s.size() == 0 || aux_f.size() == 0) {
return absl::InvalidArgumentError(absl::StrFormat(
"'%s' is '%s', which is a spline profile, so '%s_aux_s' and "
"'%s_aux_f' must be given\n",
type_key, type_name, aux_key, aux_key));
type_key, type_name, coefficient_key, coefficient_key));
}
if (aux_s.size() != aux_f.size()) {
return absl::InvalidArgumentError(absl::StrFormat(
"'%s_aux_s' and '%s_aux_f' must have the same number of entries, "
"but have %d and %d\n",
aux_key, aux_key, aux_s.size(), aux_f.size()));
coefficient_key, coefficient_key, aux_s.size(), aux_f.size()));
}
const int minimum_points = parameterization->MinimumSplinePoints();
if (aux_s.size() < minimum_points) {
return absl::InvalidArgumentError(absl::StrFormat(
"'%s' is '%s', which needs at least %d spline points, but "
"'%s_aux_s' has %d\n",
type_key, type_name, minimum_points, coefficient_key, aux_s.size()));
}
for (Eigen::Index i = 1; i < aux_s.size(); ++i) {
if (aux_s[i] <= aux_s[i - 1]) {
return absl::InvalidArgumentError(absl::StrFormat(
"'%s_aux_s' must increase strictly, but entries %d and %d are "
"%g and %g\n",
coefficient_key, i - 1, i, aux_s[i - 1], aux_s[i]));
}
}
// The Akima and cubic evaluators return zero outside their knots; the
// line segments continue their end segments instead.
const bool zero_outside_knots =
type_name.compare(0, 12, "akima_spline") == 0 ||
type_name.compare(0, 12, "cubic_spline") == 0;
if (zero_outside_knots &&
(aux_s[0] > 0.0 || aux_s[aux_s.size() - 1] < 1.0)) {
return absl::InvalidArgumentError(absl::StrFormat(
"'%s' is '%s', which is evaluated only inside its knots, so "
"'%s_aux_s' must run from 0 to 1, but runs from %g to %g\n",
type_key, type_name, coefficient_key, aux_s[0],
aux_s[aux_s.size() - 1]));
}
}

Expand Down Expand Up @@ -1425,7 +1464,7 @@ absl::Status IsConsistent(const VmecINDATA& vmec_indata,
// pmass_type, am_aux_s, am_aux_f
if (absl::Status status = CheckProfile(
"pmass_type", vmec_indata.pmass_type, ProfileType::PRESSURE, "am",
vmec_indata.am_aux_s, vmec_indata.am_aux_f);
vmec_indata.am, vmec_indata.am_aux_s, vmec_indata.am_aux_f);
!status.ok()) {
return status;
}
Expand Down Expand Up @@ -1455,17 +1494,17 @@ absl::Status IsConsistent(const VmecINDATA& vmec_indata,

// piota_type, ai_aux_s, ai_aux_f. Checked for either ncurr: piota is the
// initial guess for the iota profile even in a current-constrained run.
if (absl::Status status =
CheckProfile("piota_type", vmec_indata.piota_type, ProfileType::IOTA,
"ai", vmec_indata.ai_aux_s, vmec_indata.ai_aux_f);
if (absl::Status status = CheckProfile(
"piota_type", vmec_indata.piota_type, ProfileType::IOTA, "ai",
vmec_indata.ai, vmec_indata.ai_aux_s, vmec_indata.ai_aux_f);
!status.ok()) {
return status;
}

// pcurr_type, ac_aux_s, ac_aux_f. Ignored for ncurr == 0, still checked.
if (absl::Status status = CheckProfile(
"pcurr_type", vmec_indata.pcurr_type, ProfileType::CURRENT, "ac",
vmec_indata.ac_aux_s, vmec_indata.ac_aux_f);
vmec_indata.ac, vmec_indata.ac_aux_s, vmec_indata.ac_aux_f);
!status.ok()) {
return status;
}
Expand Down
109 changes: 109 additions & 0 deletions src/vmecpp/cpp/vmecpp/common/vmec_indata/vmec_indata_test.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#include <H5File.h>

#include <filesystem>
#include <initializer_list>
#include <sstream>
#include <string>
#include <vector>
Expand Down Expand Up @@ -267,6 +268,114 @@ TEST(TestVmecINDATA, CheckSplineProfilesNeedKnots) {
EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
}

// Below its knot count a spline evaluator returns zero for every s, so the run
// converges to an equilibrium without the requested profile. The cubic and
// Akima families need four knots, the line segments two.
TEST(TestVmecINDATA, CheckSplineProfilesNeedEnoughKnots) {
struct Case {
std::string name;
int minimum;
};
for (const Case& c : {Case{"akima_spline", 4}, Case{"cubic_spline", 4},
Case{"line_segment", 2}}) {
VmecINDATA indata;
indata.pmass_type = c.name;

indata.am_aux_s = Eigen::VectorXd::LinSpaced(c.minimum - 1, 0.0, 1.0);
indata.am_aux_f = Eigen::VectorXd::Zero(c.minimum - 1);
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok())
<< c.name;

indata.am_aux_s = Eigen::VectorXd::LinSpaced(c.minimum, 0.0, 1.0);
indata.am_aux_f = Eigen::VectorXd::Zero(c.minimum);
EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok())
<< c.name;
}

// the current and iota profiles are held to the same counts
for (const Case& c :
{Case{"akima_spline_ip", 4}, Case{"line_segment_i", 2}}) {
VmecINDATA indata;
indata.ncurr = 1;
indata.pcurr_type = c.name;
indata.ac_aux_s = Eigen::VectorXd::LinSpaced(c.minimum - 1, 0.0, 1.0);
indata.ac_aux_f = Eigen::VectorXd::Zero(c.minimum - 1);
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok())
<< c.name;
}

VmecINDATA iota_indata;
iota_indata.piota_type = "cubic_spline";
iota_indata.ai_aux_s = Eigen::VectorXd::LinSpaced(3, 0.0, 1.0);
iota_indata.ai_aux_f = Eigen::VectorXd::Zero(3);
EXPECT_FALSE(IsConsistent(iota_indata, /*enable_info_messages=*/false).ok());
}

// Below its coefficient count a closed-form evaluator returns zero for every
// s, so the run converges to an equilibrium without the requested profile.
// gauss_trunc needs 2 coefficients, two_power and two_power_gs 3 and
// two_lorentz 8; the zero-padded power series accept any count.
TEST(TestVmecINDATA, CheckClosedFormProfilesNeedEnoughCoefficients) {
struct Case {
std::string name;
int minimum;
};
for (const Case& c : {Case{"gauss_trunc", 2}, Case{"two_power", 3},
Case{"two_power_gs", 3}, Case{"two_lorentz", 8}}) {
VmecINDATA indata;
indata.pmass_type = c.name;
indata.am = Eigen::VectorXd::Ones(c.minimum - 1);
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok())
<< c.name;
indata.am = Eigen::VectorXd::Ones(c.minimum);
EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok())
<< c.name;
}

// the current profile is held to the same counts
VmecINDATA current_indata;
current_indata.ncurr = 1;
current_indata.pcurr_type = "two_power";
current_indata.ac = Eigen::VectorXd::Ones(2);
EXPECT_FALSE(
IsConsistent(current_indata, /*enable_info_messages=*/false).ok());
current_indata.ac = Eigen::VectorXd::Ones(3);
EXPECT_TRUE(
IsConsistent(current_indata, /*enable_info_messages=*/false).ok());
}

// The spline evaluators walk the knots in order and, for the Akima and cubic
// families, return zero outside them where Fortran VMEC stops; the line
// segments continue their end segments and may stop short of the radius.
TEST(TestVmecINDATA, CheckSplineKnotsIncreaseAndSpanTheRadius) {
const auto knots = [](std::initializer_list<double> values) {
Eigen::VectorXd v(static_cast<Eigen::Index>(values.size()));
Eigen::Index i = 0;
for (double value : values) v[i++] = value;
return v;
};
VmecINDATA indata;
indata.pmass_type = "cubic_spline";
indata.am_aux_f = Eigen::VectorXd::Zero(5);
indata.am_aux_s = knots({0.0, 0.25, 0.5, 0.75, 1.0});
EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
indata.am_aux_s = knots({0.0, 0.5, 0.25, 0.75, 1.0});
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
indata.am_aux_s = knots({0.0, 0.25, 0.25, 0.75, 1.0});
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
indata.am_aux_s = knots({0.0, 0.2, 0.4, 0.6, 0.9});
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
indata.am_aux_s = knots({0.1, 0.3, 0.5, 0.7, 1.0});
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok());

indata.pmass_type = "line_segment";
indata.am_aux_f = Eigen::VectorXd::Zero(3);
indata.am_aux_s = knots({0.1, 0.5, 0.9});
EXPECT_TRUE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
indata.am_aux_s = knots({0.1, 0.9, 0.5});
EXPECT_FALSE(IsConsistent(indata, /*enable_info_messages=*/false).ok());
}

TEST(TestVmecINDATA, ToJson) {
const absl::StatusOr<std::string> indata_json =
ReadFile("vmecpp/test_data/cth_like_free_bdy.json");
Expand Down
84 changes: 58 additions & 26 deletions src/vmecpp/cpp/vmecpp/vmec/profile_parameterization_data/profile_parameterization_data.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,91 +39,123 @@ std::vector<ProfileParameterizationData> BuildProfileParameterizations() {
all.reserve(NUM_PARAM);
all.emplace_back("---invalid---", /*allowedForPres=*/false,
/*allowedForCurr*/ false, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("power_series", /*allowedForPres=*/true,
/*allowedForCurr*/ true, /*allowedForIota*/ true,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("power_series_i", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("gauss_trunc", /*allowedForPres=*/true,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 2);
all.emplace_back("sum_atan", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ true,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("two_lorentz", /*allowedForPres=*/true,
/*allowedForCurr*/ false, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 8);
all.emplace_back("two_power", /*allowedForPres=*/true,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 3);
all.emplace_back("two_power_gs", /*allowedForPres=*/true,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 3);
all.emplace_back("akima_spline", /*allowedForPres=*/true,
/*allowedForCurr*/ false, /*allowedForIota*/ true,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 4,
/*minimumCoefficients*/ 0);
all.emplace_back("akima_spline_i", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 4,
/*minimumCoefficients*/ 0);
all.emplace_back("akima_spline_ip", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 4,
/*minimumCoefficients*/ 0);
all.emplace_back("cubic_spline", /*allowedForPres=*/true,
/*allowedForCurr*/ false, /*allowedForIota*/ true,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 4,
/*minimumCoefficients*/ 0);
all.emplace_back("cubic_spline_i", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 4,
/*minimumCoefficients*/ 0);
all.emplace_back("cubic_spline_ip", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 4,
/*minimumCoefficients*/ 0);
all.emplace_back("pedestal", /*allowedForPres=*/true,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("rational", /*allowedForPres=*/true,
/*allowedForCurr*/ true, /*allowedForIota*/ true,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("line_segment", /*allowedForPres=*/true,
/*allowedForCurr*/ false, /*allowedForIota*/ true,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 2,
/*minimumCoefficients*/ 0);
all.emplace_back("line_segment_i", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 2,
/*minimumCoefficients*/ 0);
all.emplace_back("line_segment_ip", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ true);
/*minimumSplinePoints*/ 2,
/*minimumCoefficients*/ 0);
all.emplace_back("nice_quadratic", /*allowedForPres=*/false,
/*allowedForCurr*/ false, /*allowedForIota*/ true,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("sum_cossq_s", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("sum_cossq_sqrts", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
all.emplace_back("sum_cossq_s_free", /*allowedForPres=*/false,
/*allowedForCurr*/ true, /*allowedForIota*/ false,
/*needsSplineData*/ false);
/*minimumSplinePoints*/ 0,
/*minimumCoefficients*/ 0);
return all;
}

} // namespace

ProfileParameterizationData::ProfileParameterizationData(
const std::string& name, bool allowedForPres, bool allowedForCurr,
bool allowedForIota, bool needsSplineData)
bool allowedForIota, int minimumSplinePoints, int minimumCoefficients)
: name_(name),
needsSplineData_(needsSplineData),
minimumSplinePoints_(minimumSplinePoints),
minimumCoefficients_(minimumCoefficients),
allowedFor_({.pres = allowedForPres,
.curr = allowedForCurr,
.iota = allowedForIota}) {}

const std::string& ProfileParameterizationData::Name() const { return name_; }

int ProfileParameterizationData::MinimumSplinePoints() const {
return minimumSplinePoints_;
}

int ProfileParameterizationData::MinimumCoefficients() const {
return minimumCoefficients_;
}

bool ProfileParameterizationData::NeedsSplineData() const {
return needsSplineData_;
return minimumSplinePoints_ > 0;
}

AllowedFor ProfileParameterizationData::IsAllowedFor() const {
Expand Down
Loading
Loading