Skip to content

[WOut] Add currumnc, currvmnc (untested) #235

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

Draft
wants to merge 1 commit into
base: graphite-base/235
Choose a base branch
from
Draft

[WOut] Add currumnc, currvmnc (untested) #235

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
9 changes: 9 additions & 0 deletions src/vmecpp/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -467,6 +467,8 @@ class VmecWOut(pydantic.BaseModel):
bsubsmns: npyd.NDArray[npyd.Shape["* dim1, * dim2"], float]
bsupumnc: npyd.NDArray[npyd.Shape["* dim1, * dim2"], float]
bsupvmnc: npyd.NDArray[npyd.Shape["* dim1, * dim2"], float]
currumnc: npyd.NDArray[npyd.Shape["* dim1, * dim2"], float]
currvmnc: npyd.NDArray[npyd.Shape["* dim1, * dim2"], float]
rmnc: npyd.NDArray[npyd.Shape["* dim1, * dim2"], float]
zmns: npyd.NDArray[npyd.Shape["* dim1, * dim2"], float]
lmns: npyd.NDArray[npyd.Shape["* mnmax, * n_surfaces"], float]
Expand Down Expand Up @@ -752,6 +754,8 @@ def save(self, out_path: str | Path) -> None:
"bsubsmns",
"bsupumnc",
"bsupvmnc",
"currumnc",
"currvmnc",
]:
fnc.createVariable(varname, np.float64, ("radius", "mn_mode_nyq"))
fnc[varname][:] = getattr(self, varname).T[:]
Expand Down Expand Up @@ -904,6 +908,8 @@ def _from_cpp_wout(cpp_wout: _vmecpp.VmecppWOut) -> VmecWOut:
attrs["rmnc"] = cpp_wout.rmnc.T
attrs["zmns"] = cpp_wout.zmns.T
attrs["bsubsmns"] = cpp_wout.bsubsmns.T
attrs["currumnc"] = cpp_wout.currumnc.T
attrs["currvmnc"] = cpp_wout.currvmnc.T

# This is a VMEC++-only quantity but it's transposed when
# stored in a wout file for consistency with lmns.
Expand Down Expand Up @@ -1046,6 +1052,9 @@ def _to_cpp_wout(self) -> _vmecpp.WOutFileContents:
cpp_wout.zmns = self.zmns.T
cpp_wout.bsubsmns = self.bsubsmns.T

cpp_wout.currumnc = self.currumnc.T
cpp_wout.currvmnc = self.currvmnc.T

# This is a VMEC++-only quantity but it's transposed when
# stored in a wout file for consistency with lmns.
cpp_wout.lmns_full = self.lmns_full.T
Expand Down
89 changes: 88 additions & 1 deletion src/vmecpp/cpp/vmecpp/vmec/output_quantities/output_quantities.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,14 @@

#include <H5Cpp.h>

#include <Eigen/Dense> // VectorXd
#include <algorithm>
#include <memory>
#include <string>
#include <vector>

#include <Eigen/Dense> // VectorXd
#include "absl/log/check.h"
#include "absl/log/log.h"
#include "util/hdf5_io/hdf5_io.h"
#include "util/testing/numerical_comparison_lib.h"
#include "vmecpp/common/util/util.h"
Expand Down Expand Up @@ -4166,6 +4167,86 @@ vmecpp::Threed1ShafranovIntegrals vmecpp::ComputeThreed1ShafranovIntegrals(
return result;
} // ComputeThreed1ShafranovIntegrals

void vmecpp::ComputeCurrents(const RowMatrixXd& bsubsmns,
const RowMatrixXd& bsubumnc,
const RowMatrixXd& bsubvmnc,
const Eigen::VectorXi& xm_nyq,
const Eigen::VectorXi& xn_nyq,
RowMatrixXd& currumnc, RowMatrixXd& currvmnc) {
const int ns = bsubsmns.rows();
Copy link

Choose a reason for hiding this comment

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

warning: narrowing conversion from 'Index' (aka 'long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]

  const int ns = bsubsmns.rows();
                 ^

const double ohs = ns - 1;
const double hs = 1.0 / ohs;
LOG(INFO) << "ComputeCurrents: ns = " << ns << ", ohs = " << ohs
<< ", hs = " << hs;
// Print dimensions of all matrices
LOG(INFO) << "bsubsmns: " << bsubsmns.rows() << " x " << bsubsmns.cols();
LOG(INFO) << "bsubumnc: " << bsubumnc.rows() << " x " << bsubumnc.cols();
LOG(INFO) << "bsubvmnc: " << bsubvmnc.rows() << " x " << bsubvmnc.cols();
LOG(INFO) << "xm_nyq: yt" << xm_nyq.size();
LOG(INFO) << "xn_nyq: " << xn_nyq.size();
LOG(INFO) << "currumnc: " << currumnc.rows() << " x " << currumnc.cols();
LOG(INFO) << "currvmnc: " << currvmnc.rows() << " x " << currvmnc.cols();
// Setup radial grid points
VectorXd shalf = VectorXd::Zero(ns);
VectorXd sfull = VectorXd::Zero(ns);
for (int js = 0; js < ns; js++) {
shalf(js) = std::sqrt(hs * (js - 1.5));
sfull(js) = std::sqrt(hs * (js - 1.0));
}

// Compute currents on interior points
for (int js = 1; js < ns - 2; js++) {
for (int mn = 0; mn < xm_nyq.size(); mn++) {
double t1{}, t2{}, t3{};
double bu0{}, bu1{}, bv0{}, bv1{};
if (static_cast<int>(xm_nyq(mn)) % 2 == 1) {
// Odd m modes
t1 = 0.5 *
(shalf(js + 1) * bsubsmns(js + 1, mn) +
shalf(js) * bsubsmns(js, mn)) /
sfull(js);

bu0 = bsubumnc(js, mn) / shalf(js);
bu1 = bsubumnc(js + 1, mn) / shalf(js + 1);
t2 = ohs * (bu1 - bu0) * sfull(js) + 0.25 * (bu0 + bu1) / sfull(js);

bv0 = bsubvmnc(js, mn) / shalf(js);
bv1 = bsubvmnc(js + 1, mn) / shalf(js + 1);
t3 = ohs * (bv1 - bv0) * sfull(js) + 0.25 * (bv0 + bv1) / sfull(js);
} else {
// Even m modes
t1 = 0.5 * (bsubsmns(js + 1, mn) + bsubsmns(js, mn));
t2 = ohs * (bsubumnc(js + 1, mn) - bsubumnc(js, mn));
t3 = ohs * (bsubvmnc(js + 1, mn) - bsubvmnc(js, mn));
}
currumnc(js, mn) = -xn_nyq(mn) * t1 - t3;
currvmnc(js, mn) = -xm_nyq(mn) * t1 + t2;
}
}
LOG(INFO) << "Completed main loop";

// Axis boundary conditions
for (int mn = 0; mn < xm_nyq.size(); mn++) {
if (xm_nyq(mn) <= 1) {
currvmnc(1, mn) = 2 * currvmnc(2, mn) - currvmnc(3, mn);
currumnc(1, mn) = 2 * currumnc(2, mn) - currumnc(3, mn);
} else {
currvmnc(1, mn) = 0;
currumnc(1, mn) = 0;
}
}

// Edge boundary conditions
for (int mn = 0; mn < xm_nyq.size(); mn++) {
currumnc(ns - 1, mn) = 2 * currumnc(ns - 2, mn) - currumnc(ns - 3, mn);
currvmnc(ns - 1, mn) = 2 * currvmnc(ns - 2, mn) - currvmnc(ns - 3, mn);
}

// Normalize by mu0
currumnc = currumnc / vmecpp::MU_0;
currvmnc = currvmnc / vmecpp::MU_0;
}

vmecpp::WOutFileContents vmecpp::ComputeWOutFileContents(
const VmecINDATA& indata, const Sizes& s, const FourierBasisFastPoloidal& t,
const FlowControl& fc, const VmecConstants& constants,
Expand Down Expand Up @@ -4609,6 +4690,12 @@ vmecpp::WOutFileContents vmecpp::ComputeWOutFileContents(
2.0 * wout.bsubsmns(1, mn_nyq) - wout.bsubsmns(2, mn_nyq);
} // mn_nyq

wout.currumnc =
RowMatrixXd::Zero(fc.ns, s.mnmax_nyq); // full-grid
wout.currvmnc =
RowMatrixXd::Zero(fc.ns, s.mnmax_nyq); // full-grid
ComputeCurrents(wout.bsubsmns, wout.bsubumnc, wout.bsubvmnc, wout.xm_nyq,
wout.xn_nyq, wout.currumnc, wout.currvmnc);
// -------------------
// non-stellarator-symmetric Fourier coefficients

Expand Down
11 changes: 11 additions & 0 deletions src/vmecpp/cpp/vmecpp/vmec/output_quantities/output_quantities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1159,6 +1159,12 @@ struct WOutFileContents {
// half-grid: contravariant B^\zeta
RowMatrixXd bsupvmnc;

// full-grid: covariant J^\theta
RowMatrixXd currumnc;

// full-grid: covariant J^\zeta
RowMatrixXd currvmnc;

// -------------------
// non-stellarator-symmetric Fourier coefficients

Expand Down Expand Up @@ -1328,6 +1334,11 @@ CovariantBDerivatives LowPassFilterCovariantB(
void ExtrapolateBSubS(const Sizes& s, const FlowControl& fc,
BSubSFull& m_bsubs_full);

void ComputeCurrents(const RowMatrixXd& bsubsmns, const RowMatrixXd& bsubumnc,
const RowMatrixXd& bsubvmnc, const Eigen::VectorXi& xm_nyq,
const Eigen::VectorXi& xn_nyq, RowMatrixXd& currumnc,
RowMatrixXd& currvmnc);

JxBOutFileContents ComputeJxBOutputFileContents(
const Sizes& s, const FlowControl& fc,
const VmecInternalResults& vmec_internal_results,
Expand Down
2 changes: 2 additions & 0 deletions 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 @@ -593,6 +593,8 @@ PYBIND11_MODULE(_vmecpp, m) {
.def_readwrite("bsubsmns_full", &vmecpp::WOutFileContents::bsubsmns_full)
.def_readwrite("bsupumnc", &vmecpp::WOutFileContents::bsupumnc)
.def_readwrite("bsupvmnc", &vmecpp::WOutFileContents::bsupvmnc)
.def_readwrite("currumnc", &vmecpp::WOutFileContents::currumnc)
.def_readwrite("currvmnc", &vmecpp::WOutFileContents::currvmnc)
//
.def_readwrite("raxis_s", &vmecpp::WOutFileContents::raxis_s)
.def_readwrite("zaxis_c", &vmecpp::WOutFileContents::zaxis_c)
Expand Down
Loading