Skip to content

Test normalize currents results in different fields #320

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
jurasic-pf wants to merge 1 commit into
base: jurasic/normalization-was-ignored
Choose a base branch
from
Open

Test normalize currents results in different fields #320

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
68 changes: 31 additions & 37 deletions src/vmecpp/cpp/vmecpp/common/makegrid_lib/makegrid_lib.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,14 +96,14 @@ absl::StatusOr<MakegridParameters> ImportMakegridParametersFromJson(

MakegridParameters makegrid_parameters;

// normalize_by_currents
// normalize_by_currents (external name for normalize_by_num_windings)
absl::StatusOr<std::optional<bool>> maybe_normalize_by_currents =
JsonReadBool(j, "normalize_by_currents");
if (!maybe_normalize_by_currents.ok()) {
return maybe_normalize_by_currents.status();
} else {
if (maybe_normalize_by_currents->has_value()) {
makegrid_parameters.normalize_by_currents =
makegrid_parameters.normalize_by_num_windings =
maybe_normalize_by_currents->value();
} else {
// MakegridParameters must be fully populated
Expand Down Expand Up @@ -444,27 +444,25 @@ absl::StatusOr<MagneticFieldResponseTable> ComputeMagneticFieldResponseTable(
sin_phi(index_phi) = std::sin(phi);
}

// Migrate num_windings into circuit currents so that normalization works
// correctly. This is done on a local copy to avoid mutating the caller's
// MagneticConfiguration.
MagneticConfiguration magnetic_config_with_currents = magnetic_configuration;
absl::Status nw_status =
NumWindingsToCircuitCurrents(magnetic_config_with_currents);
if (!nw_status.ok()) {
return nw_status;
// When normalizing, migrate num_windings into circuit currents so we can
// then set each circuit to unit current. When not normalizing, use the
// original configuration directly (ABSCAB handles num_windings natively).
MagneticConfiguration working_config = magnetic_configuration;
if (makegrid_parameters.normalize_by_num_windings) {
absl::Status nw_status = NumWindingsToCircuitCurrents(working_config);
if (!nw_status.ok()) {
return nw_status;
}
}

// Make a backup of the full vector of circuit currents (now including
// num_windings).
absl::StatusOr<std::vector<double>> maybe_original_currents =
GetCircuitCurrents(magnetic_config_with_currents);
absl::StatusOr<Eigen::VectorXd> maybe_original_currents =
GetCircuitCurrents(working_config);
if (!maybe_original_currents.ok()) {
return maybe_original_currents.status();
}
const Eigen::VectorXd& original_currents = *maybe_original_currents;

const int number_of_serial_circuits =
magnetic_config_with_currents.serial_circuits_size();
const int number_of_serial_circuits = working_config.serial_circuits_size();

MagneticFieldResponseTable response_table_b;
response_table_b.parameters = makegrid_parameters;
Expand All @@ -489,12 +487,11 @@ absl::StatusOr<MagneticFieldResponseTable> ComputeMagneticFieldResponseTable(
++circuit_index) {
// make second internal copy for being able to set inidividually only one
// circuit current to != 0
MagneticConfiguration m_magnetic_configuration =
magnetic_config_with_currents;
MagneticConfiguration m_magnetic_configuration = working_config;

Eigen::VectorXd currents_for_circuit;
currents_for_circuit.setZero(number_of_serial_circuits);
if (makegrid_parameters.normalize_by_currents) {
if (makegrid_parameters.normalize_by_num_windings) {
currents_for_circuit[circuit_index] = 1.0;
} else {
currents_for_circuit[circuit_index] = original_currents[circuit_index];
Expand Down Expand Up @@ -600,27 +597,25 @@ absl::StatusOr<MakegridCachedVectorPotential> ComputeVectorPotentialCache(
sin_phi[index_phi] = std::sin(phi);
}

// Migrate num_windings into circuit currents so that normalization works
// correctly. This is done on a local copy to avoid mutating the caller's
// MagneticConfiguration.
MagneticConfiguration magnetic_config_with_currents = magnetic_configuration;
absl::Status nw_status =
NumWindingsToCircuitCurrents(magnetic_config_with_currents);
if (!nw_status.ok()) {
return nw_status;
// When normalizing, migrate num_windings into circuit currents so we can
// then set each circuit to unit current. When not normalizing, use the
// original configuration directly (ABSCAB handles num_windings natively).
MagneticConfiguration working_config = magnetic_configuration;
if (makegrid_parameters.normalize_by_num_windings) {
absl::Status nw_status = NumWindingsToCircuitCurrents(working_config);
if (!nw_status.ok()) {
return nw_status;
}
}

// Make a backup of the full vector of circuit currents (now including
// num_windings).
absl::StatusOr<std::vector<double>> maybe_original_currents =
GetCircuitCurrents(magnetic_config_with_currents);
absl::StatusOr<Eigen::VectorXd> maybe_original_currents =
GetCircuitCurrents(working_config);
if (!maybe_original_currents.ok()) {
return maybe_original_currents.status();
}
const Eigen::VectorXd& original_currents = *maybe_original_currents;

const int number_of_serial_circuits =
magnetic_config_with_currents.serial_circuits_size();
const int number_of_serial_circuits = working_config.serial_circuits_size();

MakegridCachedVectorPotential response_table_a;
response_table_a.parameters = makegrid_parameters;
Expand All @@ -645,12 +640,11 @@ absl::StatusOr<MakegridCachedVectorPotential> ComputeVectorPotentialCache(
++circuit_index) {
// make second internal copy for being able to set inidividually only one
// circuit current to != 0
MagneticConfiguration m_magnetic_configuration =
magnetic_config_with_currents;
MagneticConfiguration m_magnetic_configuration = working_config;

Eigen::VectorXd currents_for_circuit;
currents_for_circuit.setZero(number_of_serial_circuits);
if (makegrid_parameters.normalize_by_currents) {
if (makegrid_parameters.normalize_by_num_windings) {
currents_for_circuit[circuit_index] = 1.0;
} else {
currents_for_circuit[circuit_index] = original_currents[circuit_index];
Expand Down Expand Up @@ -990,7 +984,7 @@ absl::Status WriteMakegridNetCDFFile(
CHECK_EQ(nc_put_var(ncid, id_variable_coil_group, coil_group_names.c_str()),
NC_NOERR);

if (makegrid_parameters.normalize_by_currents) {
if (makegrid_parameters.normalize_by_num_windings) {
CHECK_EQ(nc_put_var(ncid, id_variable_mgrid_mode, &kNormalizeByCurrents),
NC_NOERR);
} else {
Expand Down
2 changes: 1 addition & 1 deletion src/vmecpp/cpp/vmecpp/common/makegrid_lib/makegrid_lib.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ using RowMatrix3Xd = Eigen::Matrix<double, 3, Eigen::Dynamic, Eigen::RowMajor>;
struct MakegridParameters {
// If true, normalize the magnetic field to the coil currents and number of
// windings.
bool normalize_by_currents = false;
bool normalize_by_num_windings = false;

// If true, compute the magnetic field and vector potential only on the
// stellarator-symmetric half-period and mirror it into the other half-period.
Expand Down
79 changes: 48 additions & 31 deletions src/vmecpp/cpp/vmecpp/common/makegrid_lib/makegrid_lib_test.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,13 +68,18 @@ TEST(TestMakegridLib, CheckMakeCylindricalGridSanityChecks) {

MakegridParameters makegrid_parameters = {
// corresponding to mgrid_mode = 'R'
.normalize_by_currents = false, .assume_stellarator_symmetry = false,
.number_of_field_periods = 5, .r_grid_minimum = 1.0,
.r_grid_maximum = 2.0, .number_of_r_grid_points = 11,
.z_grid_minimum = -0.6, .z_grid_maximum = 0.6,
.number_of_z_grid_points = 13, .number_of_phi_grid_points = 18};

// both Boolean options are ok for normalize_by_currents
.normalize_by_num_windings = false,
.assume_stellarator_symmetry = false,
.number_of_field_periods = 5,
.r_grid_minimum = 1.0,
.r_grid_maximum = 2.0,
.number_of_r_grid_points = 11,
.z_grid_minimum = -0.6,
.z_grid_maximum = 0.6,
.number_of_z_grid_points = 13,
.number_of_phi_grid_points = 18};

// both Boolean options are ok for normalize_by_num_windings

// both Boolean options are ok for assume_stellarator_symmetry

Expand Down Expand Up @@ -113,17 +118,21 @@ TEST(TestMakegridLib, CheckMakeCylindricalGridSanityChecks) {
TEST(TestMakegridLib, CheckMakeCylindricalGrid) {
static constexpr double kTolerance = 1.0e-15;

MakegridParameters makegrid_parameters = {
// corresponding to mgrid_mode = 'R'
.normalize_by_currents = false, .assume_stellarator_symmetry = true,
.number_of_field_periods = 5, .r_grid_minimum = 1.0,
.r_grid_maximum = 2.0, .number_of_r_grid_points = 11,
.z_grid_minimum = -0.6, .z_grid_maximum = 0.6,
.number_of_z_grid_points = 13, .number_of_phi_grid_points = 18};
MakegridParameters makegrid_parameters = {// corresponding to mgrid_mode = 'R'
.normalize_by_num_windings = false,
.assume_stellarator_symmetry = true,
.number_of_field_periods = 5,
.r_grid_minimum = 1.0,
.r_grid_maximum = 2.0,
.number_of_r_grid_points = 11,
.z_grid_minimum = -0.6,
.z_grid_maximum = 0.6,
.number_of_z_grid_points = 13,
.number_of_phi_grid_points = 18};

// for now, make sure that the struct initialization above correctly
// identified the members
ASSERT_FALSE(makegrid_parameters.normalize_by_currents);
ASSERT_FALSE(makegrid_parameters.normalize_by_num_windings);
ASSERT_TRUE(makegrid_parameters.assume_stellarator_symmetry);
ASSERT_EQ(makegrid_parameters.number_of_field_periods, 5);
ASSERT_EQ(makegrid_parameters.r_grid_minimum, 1.0);
Expand Down Expand Up @@ -436,17 +445,21 @@ TEST(TestMakegridLib, CheckComputeMagneticFieldResponseTable) {

// NOTE: These parameters have to be consistent with the MGRID_NLI namelist
// in the `coils.test_*` input files.
MakegridParameters makegrid_parameters = {
// corresponding to mgrid_mode = 'R'
.normalize_by_currents = false, .assume_stellarator_symmetry = true,
.number_of_field_periods = 5, .r_grid_minimum = 1.0,
.r_grid_maximum = 2.0, .number_of_r_grid_points = 11,
.z_grid_minimum = -0.6, .z_grid_maximum = 0.6,
.number_of_z_grid_points = 13, .number_of_phi_grid_points = 18};
MakegridParameters makegrid_parameters = {// corresponding to mgrid_mode = 'R'
.normalize_by_num_windings = false,
.assume_stellarator_symmetry = true,
.number_of_field_periods = 5,
.r_grid_minimum = 1.0,
.r_grid_maximum = 2.0,
.number_of_r_grid_points = 11,
.z_grid_minimum = -0.6,
.z_grid_maximum = 0.6,
.number_of_z_grid_points = 13,
.number_of_phi_grid_points = 18};

// for now, make sure that the struct initialization above correctly
// identified the members
ASSERT_FALSE(makegrid_parameters.normalize_by_currents);
ASSERT_FALSE(makegrid_parameters.normalize_by_num_windings);
ASSERT_TRUE(makegrid_parameters.assume_stellarator_symmetry);
ASSERT_EQ(makegrid_parameters.number_of_field_periods, 5);
ASSERT_EQ(makegrid_parameters.r_grid_minimum, 1.0);
Expand Down Expand Up @@ -594,17 +607,21 @@ TEST(TestMakegridLib, CheckComputeVectorPotentialCache) {

// NOTE: These parameters have to be consistent with the MGRID_NLI namelist
// in the `coils.test_*` input files.
MakegridParameters makegrid_parameters = {
// corresponding to mgrid_mode = 'R'
.normalize_by_currents = false, .assume_stellarator_symmetry = true,
.number_of_field_periods = 5, .r_grid_minimum = 1.0,
.r_grid_maximum = 2.0, .number_of_r_grid_points = 11,
.z_grid_minimum = -0.6, .z_grid_maximum = 0.6,
.number_of_z_grid_points = 13, .number_of_phi_grid_points = 18};
MakegridParameters makegrid_parameters = {// corresponding to mgrid_mode = 'R'
.normalize_by_num_windings = false,
.assume_stellarator_symmetry = true,
.number_of_field_periods = 5,
.r_grid_minimum = 1.0,
.r_grid_maximum = 2.0,
.number_of_r_grid_points = 11,
.z_grid_minimum = -0.6,
.z_grid_maximum = 0.6,
.number_of_z_grid_points = 13,
.number_of_phi_grid_points = 18};

// for now, make sure that the struct initialization above correctly
// identified the members
ASSERT_FALSE(makegrid_parameters.normalize_by_currents);
ASSERT_FALSE(makegrid_parameters.normalize_by_num_windings);
ASSERT_TRUE(makegrid_parameters.assume_stellarator_symmetry);
ASSERT_EQ(makegrid_parameters.number_of_field_periods, 5);
ASSERT_EQ(makegrid_parameters.r_grid_minimum, 1.0);
Expand Down
6 changes: 2 additions & 4 deletions src/vmecpp/cpp/vmecpp/common/makegrid_lib/test_data/coils.test_symmetric_odd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,8 +68,6 @@ mirror NUL
0.827436 -0.263483 -0.364319 96
1.05079 -0.1379 -0.225085 96
1.13366 0. 0. 0 1 HF-OVF
1.266 0.0 -0.523 -16 1 HF-OVF
1.266 0.0 0.523 -16 1 HF-OVF
1.268 0.0 0.573 -36 2 TVF
1.268 0.0 -0.573 -36 2 TVF
1.266 0.0 -0.523 -16 2 HF-OVF
1.266 0.0 0.523 -16 2 HF-OVF
end
2 changes: 1 addition & 1 deletion src/vmecpp/cpp/vmecpp/free_boundary/mgrid_provider/mgrid_provider.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,7 @@ absl::Status MGridProvider::LoadFields(

nextcur = static_cast<int>(coil_currents.size());

if (mgrid_params.normalize_by_currents) {
if (mgrid_params.normalize_by_num_windings) {
mgrid_mode = "S";
} else {
mgrid_mode = "R";
Expand Down
2 changes: 1 addition & 1 deletion src/vmecpp/cpp/vmecpp/vmec/pybind11/pybind_vmec.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -672,7 +672,7 @@ PYBIND11_MODULE(_vmecpp, m) {
},
py::arg("file"))
.def_readonly("normalize_by_currents",
&makegrid::MakegridParameters::normalize_by_currents)
&makegrid::MakegridParameters::normalize_by_num_windings)
.def_readonly("assume_stellarator_symmetry",
&makegrid::MakegridParameters::assume_stellarator_symmetry)
.def_readonly("number_of_field_periods",
Expand Down
57 changes: 56 additions & 1 deletion tests/test_free_boundary.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -185,4 +185,59 @@ def test_magnetic_field_response_table_loading(makegrid_params):
)
assert response.b_p.shape == response.b_r.shape
assert response.b_z.shape == response.b_p.shape
assert response.b_r[0, 0] == pytest.approx(-9.78847973e-06, abs=1e-8, rel=1e-6)
assert response.b_r[0, 0] == pytest.approx(-1.04038999e-07, abs=1e-10, rel=1e-6)


def test_normalize_by_currents(makegrid_params: vmecpp.MakegridParameters):
makegrid_params.normalize_by_currents = True
Copy link
Copy Markdown
Collaborator

@jons-pf jons-pf May 14, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think normalize_by_currents is misleading here - I chose a bad name - since we actually aim to normalize by the number of windings -> normalize_by_num_windings?

response_normalized = vmecpp.MagneticFieldResponseTable.from_coils_file(
TEST_DATA_DIR / "coils.cth_like", makegrid_params
)

makegrid_params.normalize_by_currents = False
response_unnormalized = vmecpp.MagneticFieldResponseTable.from_coils_file(
TEST_DATA_DIR / "coils.cth_like", makegrid_params
)
assert not np.allclose(
response_normalized.b_r[0], response_unnormalized.b_r[0]
), "The two responses should be different when normalize_by_currents is True/False but were identical"
assert not np.allclose(
response_normalized.b_z[0], response_unnormalized.b_z[0]
), "The two responses should be different when normalize_by_currents is True/False but were identical"
assert not np.allclose(
response_normalized.b_p[0], response_unnormalized.b_p[0]
), "The two responses should be different when normalize_by_currents is True/False but were identical"


def test_magnetic_field_response_table_scaled_raw(makegrid_params):
assert makegrid_params.normalize_by_currents
scaled_response = vmecpp.MagneticFieldResponseTable.from_coils_file(
TEST_DATA_DIR
/ ".."
/ "common"
/ "makegrid_lib"
/ "test_data"
/ "coils.test_symmetric_odd",
makegrid_params,
)
makegrid_params.normalize_by_currents = False
raw_response = vmecpp.MagneticFieldResponseTable.from_coils_file(
TEST_DATA_DIR
/ ".."
/ "common"
/ "makegrid_lib"
/ "test_data"
/ "coils.test_symmetric_odd",
makegrid_params,
)
assert scaled_response.b_r.shape[0] == 2
assert scaled_response.b_r.shape == raw_response.b_r.shape
assert scaled_response.b_z.shape == raw_response.b_z.shape
assert scaled_response.b_p.shape == raw_response.b_p.shape
# The scaled response should be the raw response divided by the currents
np.testing.assert_allclose(
scaled_response.b_r[0] * 96, raw_response.b_r[0], atol=1e-9
)
np.testing.assert_allclose(
scaled_response.b_r[1] * -16, raw_response.b_r[1], atol=1e-9
)
Loading