Skip to content
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

Phip fix #17

Merged
merged 6 commits into from
Jul 30, 2024
Merged

Phip fix #17

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
0f2450b
Initial commit.
ejpaul Jul 9, 2024
cfb88bf
Pres added to output file.
ejpaul Jul 16, 2024
a5304a4
Added pres to output.
ejpaul Jul 16, 2024
a79fcde
Handling error in case chi not present in netcdf.
ejpaul Jul 26, 2024
87fc7b2
Merge
ejpaul Jul 30, 2024
bbfeaa5
Matt comments.
ejpaul Jul 30, 2024
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
19 changes: 15 additions & 4 deletions src/_booz_xform/bindings.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,8 @@ need to call this function directly.
For stellarator-symmetric configurations this array is ignored and need not
be specified.
:param phip: Phip on vmec's full grid. (Defaults to empty Vector)
:param chi: chi on vmec's full grid. (Defaults to empty Vector)
:param pres: pressure on vmec's full grid. (Defaults to empty Vector)
)",
"ns"_a,
"iotas"_a,
Expand All @@ -76,7 +78,9 @@ need to call this function directly.
"bsubumns0"_a,
"bsubvmnc0"_a,
"bsubvmns0"_a,
"phips"_a=defaultInitPtr)
"phips"_a=defaultInitPtr,
"chi"_a=defaultInitPtr,
"pres"_a=defaultInitPtr)

.def("run", &Booz_xform::run, R"(
Run the transformation to Boozer coordinates on all the selected
Expand Down Expand Up @@ -177,10 +181,17 @@ input radial surfaces.)")
(1D float array of length ns_in, input) Toroidal flux normalized by 2*pi,
evaluated on all the magnetic surfaces for which input data was provided.)")

.def_readwrite("phip", &Booz_xform::phip, R"(
.def_readwrite("chi", &Booz_xform::chi, R"(
(1D float array of length ns_in, input) Poloidal flux normalized by 2*pi,
evaluated on all the magnetic surfaces for which input data was provided.)")

.def_readwrite("phip", &Booz_xform::phip, R"(
(1D float array of length ns_in, input) Toroidal flux normalized by 2*pi,
evaluated on all the magnetic surfaces for which input data was provided.)")

.def_readwrite("pres", &Booz_xform::pres, R"(
(1D float array of length ns_in, input) Pressure on full vmec grid.)")

.def_readwrite("rmnc", &Booz_xform::rmnc, R"(
(2D float array of size mnmax x ns_in, input) cos(m * theta_0 - n *
zeta_0) Fourier modes of the major radius R.)")
Expand Down Expand Up @@ -311,12 +322,12 @@ quantity is zero so the array will have size 0 x 0.)")
.def_readonly("gmnc_b", &Booz_xform::gmnc_b, R"(
(2D float array of size mnboz x ns_b, output) cos(m * theta_B - n *
zeta_B) Fourier modes (with respect to Boozer coordinates) of the
Jacobian of (psi, theta_B, zeta_B) coordinates.)")
Jacobian of (s, theta_B, zeta_B) coordinates.)")
ejpaul marked this conversation as resolved.
Show resolved Hide resolved

.def_readonly("gmns_b", &Booz_xform::gmns_b, R"(
(2D float array of size mnboz x ns_b, output) sin(m * theta_B - n *
zeta_B) Fourier modes (with respect to Boozer coordinates) of the
Jacobian of (psi, theta_B, zeta_B) coordinates. If the configuration is
Jacobian of (s, theta_B, zeta_B) coordinates. If the configuration is
ejpaul marked this conversation as resolved.
Show resolved Hide resolved
stellarator-symmetric, this quantity is zero so this array will have
size 0 x 0.)")

Expand Down
18 changes: 14 additions & 4 deletions src/_booz_xform/booz_xform.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,13 +304,13 @@ namespace booz_xform {

/** (size mnboz x ns_b, output) cos(m * theta_B - n * zeta_B)
Fourier modes (with respect to Boozer coordinates) of the
Jacobian of (psi, theta_B, zeta_B) coordinates.
Jacobian of (s, theta_B, zeta_B) coordinates.
ejpaul marked this conversation as resolved.
Show resolved Hide resolved
*/
Matrix gmnc_b;

/** (size mnboz x ns_b, output) sin(m * theta_B - n * zeta_B)
Fourier modes (with respect to Boozer coordinates) of the
Jacobian of (psi, theta_B, zeta_B) coordinates. If the
Jacobian of (s, theta_B, zeta_B) coordinates. If the
ejpaul marked this conversation as resolved.
Show resolved Hide resolved
configuration is stellarator-symmetric, this quantity is zero
so this array will have size 0 x 0.
*/
Expand Down Expand Up @@ -391,10 +391,18 @@ namespace booz_xform {
*/
Vector phi;

/** (size ns_in + 1, output) The derivative of the toroidal flux (not divided by (2*pi))
with respect to s, evaluated on full vmec grid.
*/
Vector phip;

/** (size ns_in + 1, output) Uniformly spaced grid going from 0 to the boundary
poloidal flux (not divided by (2*pi)), evaluated on full vmec grid.
*/
Vector phip;
Vector chi;

/** Pressure evaluated on the full vmec grid **/
Vector pres;

//! Constructor
/**
Expand Down Expand Up @@ -434,7 +442,9 @@ namespace booz_xform {
Matrix& bsubumns,
Matrix& bsubvmnc,
Matrix& bsubvmns,
Vector& phips=defaultInitPtr);
Vector& phips=defaultInitPtr,
Vector& chi=defaultInitPtr,
Vector& pres=defaultInitPtr);

//! Carry out the transformation calculation
/**
Expand Down
14 changes: 12 additions & 2 deletions src/_booz_xform/init_from_vmec.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,9 @@ void Booz_xform::init_from_vmec(int ns,
Matrix& bsubumns0,
Matrix& bsubvmnc0,
Matrix& bsubvmns0,
Vector& phip0) {
Vector& phip0,
Vector& chi0,
Vector& pres0) {
int j, k;
ns_in = ns - 1;

Expand All @@ -34,6 +36,8 @@ void Booz_xform::init_from_vmec(int ns,
bool skip_phip0 = (phip0.size()==0);
if (!skip_phip0) {
if (phip0.size() != ns) throw std::runtime_error("phip0.size() is not ns");
if (chi0.size() != ns) throw std::runtime_error("chi0.size() is not ns");
if (pres0.size() != ns) throw std::runtime_error("pres0.size() is not ns");
}
if (xm.size() != mnmax) throw std::runtime_error("Size of xm is not mnmax");
if (xn.size() != mnmax) throw std::runtime_error("Size of xn is not mnmax");
Expand Down Expand Up @@ -94,11 +98,17 @@ void Booz_xform::init_from_vmec(int ns,
}
if (!skip_phip0) {
phip.resize(ns_in + 1);
chi.resize(ns_in + 1);
pres.resize(ns_in + 1);
for (j = 0; j < ns_in + 1; j++) {
phip[j] = phip0[j];
phip[j] = -phip0[j]/(2*pi);
chi[j] = chi0[j];
pres[j] = pres0[j];
}
} else {
phip.resize(0);
chi.resize(0);
pres.resize(0);
}
// By default, prepare to do the Boozer transformation at all
// half-grid surfaces:
Expand Down
7 changes: 7 additions & 0 deletions src/_booz_xform/read_boozmn.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,11 +40,18 @@ void Booz_xform::read_boozmn(std::string filename) {
Boozer_I_all.resize(ns_in);
phi.resize(ns_in+1);
phip.resize(ns_in+1);
chi.resize(ns_in+1);
pres.resize(ns_in+1);
nc.get("iota_b", iota_in);
nc.get("bvco_b", Boozer_G_in);
nc.get("buco_b", Boozer_I_in);
nc.get("phi_b", phi);
nc.get("phip_b", phip);
try {
nc.get("chi_b", chi);
} catch (std::runtime_error error) {
}
nc.get("pres_b", pres);
for (j = 0; j < ns_in; j++) {
iota[j] = iota_in[j+1];
Boozer_G_all[j] = Boozer_G_in[j+1];
Expand Down
10 changes: 8 additions & 2 deletions src/_booz_xform/read_wout.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,9 +28,15 @@ void Booz_xform::read_wout(std::string filename, bool flux) {
nc.get("aspect", aspect);

Vector phip0;
Vector chi0;
Vector pres0;
if (flux) {
phip0.resize(ns);
nc.get("chi", phip0);
chi0.resize(ns);
pres0.resize(ns);
nc.get("chi", chi0);
nc.get("phipf", phip0);
nc.get("pres",pres0);
}

phi.resize(ns);
Expand Down Expand Up @@ -100,7 +106,7 @@ void Booz_xform::read_wout(std::string filename, bool flux) {
if (flux) {
init_from_vmec(ns, iotas, rmnc0, rmns0, zmnc0, zmns0,
lmnc0, lmns0, bmnc0, bmns0,
bsubumnc0, bsubumns0, bsubvmnc0, bsubvmns0, phip0);
bsubumnc0, bsubumns0, bsubvmnc0, bsubvmns0, phip0, chi0, pres0);
} else {
init_from_vmec(ns, iotas, rmnc0, rmns0, zmnc0, zmns0,
lmnc0, lmns0, bmnc0, bmns0,
Expand Down
13 changes: 9 additions & 4 deletions src/_booz_xform/write_boozmn.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,16 +81,21 @@ void Booz_xform::write_boozmn(std::string filename) {
// available here. Just write vectors of zeros for
// backwards-compatibility of the boozmn.nc files.
Vector pres_b(ns_in + 1), beta_b(ns_in + 1);
pres_b.setZero();
nc.put(radius_dim, "pres_b", pres_b, placeholder_str, "dimensionless");
beta_b.setZero();
nc.put(radius_dim, "beta_b", pres_b, placeholder_str, "dimensionless");
Vector phip_dummy(ns_in + 1);
Vector phip_dummy(ns_in + 1), chi_dummy(ns_in + 1);
phip_dummy.setZero();
chi_dummy.setZero();
if (phip.size()==0) {
pres_b.setZero();
nc.put(radius_dim, "phip_b", phip_dummy, placeholder_str, "Tesla * meter^2");
nc.put(radius_dim, "chi_b", chi_dummy, placeholder_str, "Tesla * meter^2");
nc.put(radius_dim, "pres_b", pres_b, placeholder_str, "dimensionless");
} else {
nc.put(radius_dim, "phip_b", phip, "Uniformly spaced grid going from 0 to the boundary poloidal flux (not divided by (2*pi)). This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");
phip[0] = 0;
nc.put(radius_dim, "phip_b", phip, "Derivative of toroidal flux (not divided by (2*pi)) with respect to s. This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");
nc.put(radius_dim, "chi_b", chi, "Uniformly spaced grid going from 0 to the boundary poloidal flux (not divided by (2*pi)). This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");
nc.put(radius_dim, "pres_b", pres, "Pressure on full vmec grid. This grid generally does not correspond to the radial grid used for other quantities in this file!", "Pascal");
}
nc.put(radius_dim, "phi_b", phi, "Uniformly spaced grid going from 0 to the boundary toroidal flux (not divided by (2*pi)). This grid generally does not correspond to the radial grid used for other quantities in this file!", "Tesla * meter^2");

Expand Down
107 changes: 56 additions & 51 deletions tests/test_regression.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,65 +22,70 @@ def test_regression(self):
boozmn_new_filename = 'boozmn_new_' + configuration + '.nc'
f = netcdf_file(os.path.join(TEST_DIR, boozmn_filename),
'r', mmap=False)
b = Booz_xform()
b.read_wout(os.path.join(TEST_DIR, wout_filename))
# Transfer parameters from the reference file to the new
# calculation
b.mboz = f.variables['mboz_b'][()]
b.nboz = f.variables['nboz_b'][()]
b.compute_surfs = f.variables['jlist'][()] - 2

b.run()
for flux in [True,False]:
b = Booz_xform()
b.read_wout(os.path.join(TEST_DIR, wout_filename),flux=flux)
# Transfer parameters from the reference file to the new
# calculation
b.mboz = f.variables['mboz_b'][()]
b.nboz = f.variables['nboz_b'][()]
b.compute_surfs = f.variables['jlist'][()] - 2

# Compare 2D arrays
vars = ['bmnc_b', 'rmnc_b', 'zmns_b', 'numns_b', 'gmnc_b']
asym = bool(f.variables['lasym__logical__'][()])
if asym:
vars += ['bmns_b', 'rmns_b', 'zmnc_b', 'numnc_b', 'gmns_b']
b.run()

rtol = 1e-12
atol = 1e-12
for var in vars:
# gmnc_b is misspelled in the fortran version
var_ref = var
if var == 'gmnc_b':
var_ref = 'gmn_b'
# Compare 2D arrays
vars = ['bmnc_b', 'rmnc_b', 'zmns_b', 'numns_b', 'gmnc_b']
asym = bool(f.variables['lasym__logical__'][()])
if asym:
vars += ['bmns_b', 'rmns_b', 'zmnc_b', 'numnc_b', 'gmns_b']

# Handle the issue that we now use the variable nu,
# whereas the boozmn format uses the variable
# p = -nu.
sign = 1
if var[:2] == 'nu':
sign = -1
var_ref = 'p' + var[2:]
rtol = 1e-12
atol = 1e-12
for var in vars:
# gmnc_b is misspelled in the fortran version
var_ref = var
if var == 'gmnc_b':
var_ref = 'gmn_b'

# Reference values:
arr1 = f.variables[var_ref][()]
# Newly computed values:
arr2 = getattr(b, var).transpose()
# Handle the issue that we now use the variable nu,
# whereas the boozmn format uses the variable
# p = -nu.
sign = 1
if var[:2] == 'nu':
sign = -1
var_ref = 'p' + var[2:]

print('abs diff in ' + var + ':', np.max(np.abs(arr1 - sign * arr2)))
np.testing.assert_allclose(arr1, sign * arr2,
rtol=rtol, atol=atol)
# Reference values:
arr1 = f.variables[var_ref][()]
# Newly computed values:
arr2 = getattr(b, var).transpose()

# Now compare some values written to the boozmn files.
b.write_boozmn(boozmn_new_filename)
f2 = netcdf_file(boozmn_new_filename)
print('abs diff in ' + var + ':', np.max(np.abs(arr1 - sign * arr2)))
np.testing.assert_allclose(arr1, sign * arr2,
rtol=rtol, atol=atol)

vars = f.variables.keys()
# These variables will not match:
exclude = ['rmax_b', 'rmin_b', 'betaxis_b', 'version', 'pres_b', 'beta_b', 'phip_b']
for var in vars:
if var in exclude:
continue
# Reference values:
arr1 = f.variables[var][()]
# Newly computed values:
arr2 = f2.variables[var][()]
print('abs diff in ' + var + ':', np.max(np.abs(arr1 - arr2)))
np.testing.assert_allclose(arr1, arr2,
rtol=rtol, atol=atol)

# Now compare some values written to the boozmn files.
b.write_boozmn(boozmn_new_filename)
f2 = netcdf_file(boozmn_new_filename)

vars = f.variables.keys()
# These variables will not match:
if flux:
exclude = ['rmax_b', 'rmin_b', 'betaxis_b', 'version', 'beta_b']
else:
exclude = ['rmax_b', 'rmin_b', 'betaxis_b', 'version', 'pres_b', 'beta_b', 'phip_b']
for var in vars:
if var in exclude:
continue
# Reference values:
arr1 = f.variables[var][()]
# Newly computed values:
arr2 = f2.variables[var][()]
print('abs diff in ' + var + ':', np.max(np.abs(arr1 - arr2)))
np.testing.assert_allclose(arr1, arr2,
rtol=rtol, atol=atol)

f.close()

if __name__ == '__main__':
Expand Down
4 changes: 4 additions & 0 deletions tests/test_write_read.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,9 +63,13 @@ def test_write_read(self):
np.testing.assert_allclose(b1.phi,b2.phi, rtol=rtol, atol=atol)
if flux:
np.testing.assert_allclose(b1.phip,b2.phip, rtol=rtol, atol=atol)
np.testing.assert_allclose(b1.chi,b2.chi, rtol=rtol, atol=atol)
np.testing.assert_allclose(b1.pres,b2.pres, rtol=rtol, atol=atol)
else:
np.testing.assert_equal(b2.phip,0)
assert len(b1.phip) == 0
np.testing.assert_equal(b2.chi,0)
assert len(b1.chi) == 0
os.remove(boozmn_filename)

if __name__ == '__main__':
Expand Down
Loading