Skip to content

Commit ca6ddfc

Browse files
jonathanschillingjurasic-pf
jonathanschilling
authored and
jurasic-pf
committed
add free-boundary method that only uses magnetic field from external coils for pure vacuum case
1 parent 40393aa commit ca6ddfc

14 files changed

Lines changed: 955 additions & 9 deletions

File tree

‎examples/data/w7x_free_bdy_vac.json‎

Lines changed: 712 additions & 0 deletions
Large diffs are not rendered by default.

‎src/vmecpp/__init__.py‎

Lines changed: 50 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,42 @@ class RestartReason(enum.Enum):
7474
not overlapping yet)"""
7575

7676

77+
class FreeBoundaryMethod(str, enum.Enum):
78+
"""Method for handling free-boundary conditions."""
79+
80+
NESTOR = "nestor"
81+
"""NEumann Solver for TORoidal systems."""
82+
83+
ONLY_COILS = "only_coils"
84+
"""Use just the coils field for the free-boundary force contribution.
85+
86+
This can be particularly useful for verification calculations.
87+
88+
Warning: This is only valid for vacuum calculations and will ignore
89+
the plasma current contribution!
90+
"""
91+
92+
BIEST = "biest"
93+
"""Boundary Integral Equation Solver for Toroidal systems."""
94+
95+
96+
def _validate_free_boundary_method(
97+
value: _vmecpp.FreeBoundaryMethod | str | FreeBoundaryMethod,
98+
) -> FreeBoundaryMethod:
99+
"""Convert various representations to FreeBoundaryMethod."""
100+
if isinstance(value, FreeBoundaryMethod):
101+
return value
102+
if isinstance(value, _vmecpp.FreeBoundaryMethod):
103+
return FreeBoundaryMethod[value.name] # type: ignore
104+
if isinstance(value, str):
105+
return FreeBoundaryMethod(value)
106+
msg = (
107+
f"Cannot convert {type(value).__name__} to FreeBoundaryMethod. "
108+
"Expected a string or FreeBoundaryMethod instance."
109+
)
110+
raise ValueError(msg)
111+
112+
77113
# This is a pure Python equivalent of VmecINDATAPyWrapper.
78114
# In the future VmecINDATAPyWrapper and the C++ VmecINDATA will merge into one type,
79115
# and this will become a Python wrapper around the one C++ VmecINDATA type.
@@ -262,6 +298,13 @@ class VmecInput(BaseModelWithNumpy):
262298
nvacskip: int = 1
263299
"""Number of iterations between full vacuum calculations."""
264300

301+
free_boundary_method: typing.Annotated[
302+
FreeBoundaryMethod,
303+
pydantic.BeforeValidator(_validate_free_boundary_method),
304+
pydantic.Field(),
305+
] = FreeBoundaryMethod.NESTOR
306+
"""Method for handling free-boundary conditions."""
307+
265308
nstep: int = 10
266309
"""Printout interval at which convergence progress is logged."""
267310

@@ -506,10 +549,15 @@ def _to_cpp_vmecindata(self) -> _vmecpp.VmecINDATA:
506549
}
507550

508551
for attr in VmecInput.model_fields:
509-
if attr in readonly_attrs:
552+
if attr in readonly_attrs or attr == "free_boundary_method":
510553
continue # these must be set separately
511554
setattr(cpp_indata, attr, getattr(self, attr))
512555

556+
# Convert Python enum to C++ enum
557+
cpp_indata.free_boundary_method = getattr(
558+
_vmecpp.FreeBoundaryMethod, self.free_boundary_method.name
559+
)
560+
513561
# this also resizes the readonly_attrs
514562
cpp_indata._set_mpol_ntor(self.mpol, self.ntor)
515563
for attr in readonly_attrs - {"mpol", "ntor"}:
@@ -2119,5 +2167,6 @@ def populate_raw_profile(
21192167
"Threed1Volumetrics",
21202168
"MakegridParameters",
21212169
"MagneticFieldResponseTable",
2170+
"FreeBoundaryMethod",
21222171
"populate_raw_profile",
21232172
]

‎src/vmecpp/cpp/vmecpp/common/vmec_indata/vmec_indata.cc‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -51,6 +51,8 @@ absl::StatusOr<FreeBoundaryMethod> FreeBoundaryMethodFromString(
5151
const std::string& free_boundary_method_string) {
5252
if (free_boundary_method_string == "nestor") {
5353
return FreeBoundaryMethod::NESTOR;
54+
} else if (free_boundary_method_string == "only_coils") {
55+
return FreeBoundaryMethod::ONLY_COILS;
5456
} else if (free_boundary_method_string == "biest") {
5557
return FreeBoundaryMethod::BIEST;
5658
}
@@ -63,6 +65,8 @@ std::string ToString(FreeBoundaryMethod free_boundary_method) {
6365
switch (free_boundary_method) {
6466
case FreeBoundaryMethod::NESTOR:
6567
return "nestor";
68+
case FreeBoundaryMethod::ONLY_COILS:
69+
return "only_coils";
6670
case FreeBoundaryMethod::BIEST:
6771
return "biest";
6872
default:
@@ -1388,10 +1392,12 @@ absl::Status IsConsistent(const VmecINDATA& vmec_indata,
13881392
// For the current state of the code, we only accept NESTOR,
13891393
// but in the future [TODO(jons)] also all other (implemented) enum values
13901394
// are valid.
1391-
if (vmec_indata.free_boundary_method != FreeBoundaryMethod::NESTOR) {
1392-
return absl::InvalidArgumentError(absl::StrFormat(
1393-
"input variable 'free_boundary_method' must be 'nestor', but is %s\n",
1394-
ToString(vmec_indata.free_boundary_method)));
1395+
if (vmec_indata.free_boundary_method != FreeBoundaryMethod::NESTOR &&
1396+
vmec_indata.free_boundary_method != FreeBoundaryMethod::ONLY_COILS) {
1397+
return absl::InvalidArgumentError(
1398+
absl::StrFormat("input variable 'free_boundary_method' must be "
1399+
"'nestor' or 'only_coils', but is %s\n",
1400+
ToString(vmec_indata.free_boundary_method)));
13951401
}
13961402
}
13971403

‎src/vmecpp/cpp/vmecpp/common/vmec_indata/vmec_indata.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,10 @@ enum class FreeBoundaryMethod : std::uint8_t {
3030
// for the free-boundary force contribution
3131
NESTOR,
3232

33+
// use just the coils field
34+
// for the free-boundary force contribution
35+
ONLY_COILS,
36+
3337
// use the Boundary Integral Equation Solver for Toroidal systems
3438
// for the free-boundary force contribution
3539
BIEST

‎src/vmecpp/cpp/vmecpp/free_boundary/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ add_subdirectory(free_boundary_base)
33
add_subdirectory(laplace_solver)
44
add_subdirectory(mgrid_provider)
55
add_subdirectory(nestor)
6+
add_subdirectory(only_coils)
67
add_subdirectory(regularized_integrals)
78
add_subdirectory(singular_integrals)
89
add_subdirectory(surface_geometry)

‎src/vmecpp/cpp/vmecpp/free_boundary/only_coils/BUILD.bazel‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
# SPDX-FileCopyrightText: 2024-present Proxima Fusion GmbH <info@proximafusion.com>
2+
#
3+
# SPDX-License-Identifier: MIT
4+
cc_library(
5+
name = "only_coils",
6+
srcs = ["only_coils.cc"],
7+
hdrs = ["only_coils.h"],
8+
visibility = ["//visibility:public"],
9+
deps = [
10+
"//vmecpp/common/util:util",
11+
"//vmecpp/common/sizes:sizes",
12+
"//vmecpp/common/fourier_basis_fast_toroidal",
13+
"//vmecpp/free_boundary/tangential_partitioning:tangential_partitioning",
14+
"//vmecpp/free_boundary/external_magnetic_field:external_magnetic_field",
15+
"//vmecpp/free_boundary/singular_integrals:singular_integrals",
16+
"//vmecpp/free_boundary/regularized_integrals:regularized_integrals",
17+
"//vmecpp/free_boundary/laplace_solver:laplace_solver",
18+
"//vmecpp/free_boundary/free_boundary_base:free_boundary_base",
19+
],
20+
)

‎src/vmecpp/cpp/vmecpp/free_boundary/only_coils/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
list (APPEND vmecpp_sources
2+
${CMAKE_CURRENT_SOURCE_DIR}/only_coils.cc
3+
${CMAKE_CURRENT_SOURCE_DIR}/only_coils.h
4+
)
5+
set (vmecpp_sources "${vmecpp_sources}" PARENT_SCOPE)

‎src/vmecpp/cpp/vmecpp/free_boundary/only_coils/only_coils.cc‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
// SPDX-FileCopyrightText: 2024-present Proxima Fusion GmbH
2+
// <info@proximafusion.com>
3+
//
4+
// SPDX-License-Identifier: MIT
5+
#include "vmecpp/free_boundary/only_coils/only_coils.h"
6+
7+
namespace vmecpp {
8+
9+
OnlyCoils::OnlyCoils(const Sizes* s, const TangentialPartitioning* tp,
10+
const MGridProvider* mgrid, std::span<double> bSqVacShare,
11+
std::span<double> vacuum_b_r_share,
12+
std::span<double> vacuum_b_phi_share,
13+
std::span<double> vacuum_b_z_share)
14+
: FreeBoundaryBase(s, tp, mgrid, bSqVacShare, vacuum_b_r_share,
15+
vacuum_b_phi_share, vacuum_b_z_share) {} // OnlyCoils
16+
17+
bool OnlyCoils::update(
18+
const std::span<const double> rCC, const std::span<const double> rSS,
19+
const std::span<const double> rSC, const std::span<const double> rCS,
20+
const std::span<const double> zSC, const std::span<const double> zCS,
21+
const std::span<const double> zCC, const std::span<const double> zSS,
22+
int signOfJacobian, const std::span<const double> rAxis,
23+
const std::span<const double> zAxis, double* bSubUVac, double* bSubVVac,
24+
double netToroidalCurrent, int ivacskip,
25+
const VmecCheckpoint& vmec_checkpoint, bool at_checkpoint_iteration) {
26+
// only need surface geometry, not all derived quantities
27+
bool full_update = false;
28+
29+
sg_.update(rCC, rSS, rSC, rCS, zSC, zCS, zCC, zSS, signOfJacobian,
30+
full_update);
31+
32+
// blindly assume netToroidalCurrent == 0.0,
33+
// since checked for that during initialization
34+
ef_.update(rAxis, zAxis, 0.0);
35+
36+
// compute net covariant magnetic field components on surface
37+
double local_bsubuvac = 0.0;
38+
double local_bsubvvac = 0.0;
39+
for (int kl = tp_.ztMin; kl < tp_.ztMax; ++kl) {
40+
int l = kl / s_.nZeta;
41+
local_bsubuvac += ef_.bSubU[kl - tp_.ztMin] * s_.wInt[l];
42+
local_bsubvvac += ef_.bSubV[kl - tp_.ztMin] * s_.wInt[l];
43+
}
44+
local_bsubuvac *= signOfJacobian * 2.0 * M_PI;
45+
46+
#ifdef _OPENMP
47+
#pragma omp single
48+
#endif // _OPENMP
49+
{
50+
*bSubUVac = 0.0;
51+
*bSubVVac = 0.0;
52+
}
53+
#ifdef _OPENMP
54+
#pragma omp barrier
55+
#endif // _OPENMP
56+
57+
#ifdef _OPENMP
58+
#pragma omp critical
59+
#endif // _OPENMP
60+
{
61+
*bSubUVac += local_bsubuvac;
62+
*bSubVVac += local_bsubvvac;
63+
}
64+
#ifdef _OPENMP
65+
#pragma omp barrier
66+
#endif // _OPENMP
67+
68+
// compute magnetic pressure from only coils: |B|^2/2
69+
for (int kl = tp_.ztMin; kl < tp_.ztMax; ++kl) {
70+
// cylindrical components of vacuum magnetic field
71+
vacuum_b_r_share_[kl] = ef_.interpBr[kl - tp_.ztMin];
72+
vacuum_b_phi_share_[kl] = ef_.interpBp[kl - tp_.ztMin];
73+
vacuum_b_z_share_[kl] = ef_.interpBz[kl - tp_.ztMin];
74+
75+
// magnetic pressure from vacuum: |B|^2/2
76+
bSqVacShare[kl] = 0.5 * (vacuum_b_r_share_[kl] * vacuum_b_r_share_[kl] +
77+
vacuum_b_phi_share_[kl] * vacuum_b_phi_share_[kl] +
78+
vacuum_b_z_share_[kl] * vacuum_b_z_share_[kl]);
79+
} // kl
80+
81+
// ... done ...
82+
83+
#ifdef _OPENMP
84+
#pragma omp barrier
85+
#endif // _OPENMP
86+
87+
// TODO(jons): could move bSubUVac, bSubVVac collection here to spare on
88+
// barrier
89+
90+
return false;
91+
} // update
92+
93+
} // namespace vmecpp

‎src/vmecpp/cpp/vmecpp/free_boundary/only_coils/only_coils.h‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
// SPDX-FileCopyrightText: 2024-present Proxima Fusion GmbH
2+
// <info@proximafusion.com>
3+
//
4+
// SPDX-License-Identifier: MIT
5+
#ifndef VMECPP_FREE_BOUNDARY_ONLY_COILS_ONLY_COILS_H_
6+
#define VMECPP_FREE_BOUNDARY_ONLY_COILS_ONLY_COILS_H_
7+
8+
#include <span>
9+
10+
#include "vmecpp/common/sizes/sizes.h"
11+
#include "vmecpp/common/util/util.h"
12+
#include "vmecpp/free_boundary/free_boundary_base/free_boundary_base.h"
13+
#include "vmecpp/free_boundary/mgrid_provider/mgrid_provider.h"
14+
#include "vmecpp/free_boundary/tangential_partitioning/tangential_partitioning.h"
15+
16+
namespace vmecpp {
17+
18+
class OnlyCoils : public FreeBoundaryBase {
19+
public:
20+
OnlyCoils(const Sizes* s, const TangentialPartitioning* tp,
21+
const MGridProvider* mgrid, std::span<double> bSqVacShare,
22+
std::span<double> vacuum_b_r_share,
23+
std::span<double> vacuum_b_phi_share,
24+
std::span<double> vacuum_b_z_share);
25+
26+
bool update(
27+
const std::span<const double> rCC, const std::span<const double> rSS,
28+
const std::span<const double> rSC, const std::span<const double> rCS,
29+
const std::span<const double> zSC, const std::span<const double> zCS,
30+
const std::span<const double> zCC, const std::span<const double> zSS,
31+
int signOfJacobian, const std::span<const double> rAxis,
32+
const std::span<const double> zAxis, double* bSubUVac, double* bSubVVac,
33+
double netToroidalCurrent, int ivacskip,
34+
const VmecCheckpoint& vmec_checkpoint = VmecCheckpoint::NONE,
35+
bool at_checkpoint_iteration = false) final;
36+
};
37+
38+
} // namespace vmecpp
39+
40+
#endif // VMECPP_FREE_BOUNDARY_ONLY_COILS_ONLY_COILS_H_

‎src/vmecpp/cpp/vmecpp/vmec/pybind11/pybind_vmec.cc‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -196,6 +196,7 @@ PYBIND11_MODULE(_vmecpp, m) {
196196

197197
py::enum_<vmecpp::FreeBoundaryMethod>(m, "FreeBoundaryMethod")
198198
.value("NESTOR", vmecpp::FreeBoundaryMethod::NESTOR)
199+
.value("ONLY_COILS", vmecpp::FreeBoundaryMethod::ONLY_COILS)
199200
.value("BIEST", vmecpp::FreeBoundaryMethod::BIEST);
200201

201202
py::enum_<vmecpp::IterationStyle>(m, "IterationStyle")

0 commit comments

Comments
 (0)