Skip to content

Commit

Permalink
Move implementation of TokamakCoordinates class out of header file
Browse files Browse the repository at this point in the history
  • Loading branch information
tomchapman committed Dec 16, 2024
1 parent 1e2a981 commit 01e5816
Show file tree
Hide file tree
Showing 3 changed files with 77 additions and 61 deletions.
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -361,6 +361,7 @@ set(BOUT_SOURCES
${CMAKE_CURRENT_BINARY_DIR}/include/bout/revision.hxx
${CMAKE_CURRENT_BINARY_DIR}/include/bout/version.hxx
./include/bout/tokamak_coordinates.hxx
./src/mesh/tokamak_coordinates.cxx
)


Expand Down
68 changes: 7 additions & 61 deletions include/bout/tokamak_coordinates.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,11 @@
#include "bout/utils.hxx"
#include "bout/vector2d.hxx"


class TokamakCoordinates {

private:

Mesh& mesh_m;
Field2D Rxy_m;
Field2D Bpxy_m;
Expand All @@ -22,73 +24,17 @@ private:
FieldMetric ShearFactor_m;
FieldMetric dx_m;

void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) {
void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor);

Rxy_m /= Lbar;
Bpxy_m /= Bbar;
Btxy_m /= Bbar;
Bxy_m /= Bbar;
hthe_m /= Lbar;
ShearFactor_m *= Lbar * Lbar * Bbar * ShearFactor;
dx_m /= Lbar * Lbar * Bbar;
}
BoutReal get_sign_of_bp();

public:
TokamakCoordinates(Mesh& mesh) : mesh_m(mesh) {

mesh.get(Rxy_m, "Rxy");
// mesh->get(Zxy, "Zxy");
mesh.get(Bpxy_m, "Bpxy");
mesh.get(Btxy_m, "Btxy");
mesh.get(Bxy_m, "Bxy");
mesh.get(hthe_m, "hthe");
mesh.get(ShearFactor_m, "sinty");
mesh.get(dx_m, "dpsi");
}
public:

BoutReal get_sign_of_bp() {
if (min(Bpxy_m, true) < 0.0) {
return -1.0;
}
return 1.0;
}
TokamakCoordinates(Mesh& mesh);

Coordinates* make_coordinates(const bool noshear, BoutReal Lbar, BoutReal Bbar,
BoutReal ShearFactor = 1.0) {

normalise(Lbar, Bbar, ShearFactor);

const BoutReal sign_of_bp = get_sign_of_bp();

if (noshear) {
ShearFactor_m = 0.0;
}

auto* coord = mesh_m.getCoordinates();

const auto g11 = SQ(Rxy_m * Bpxy_m);
const auto g22 = 1.0 / SQ(hthe_m);
const auto g33 = SQ(ShearFactor_m) * g11 + SQ(Bxy_m) / g11;
const auto g12 = 0.0;
const auto g13 = -ShearFactor_m * g11;
const auto g23 = -sign_of_bp * Btxy_m / (hthe_m * Bpxy_m * Rxy_m);

const auto g_11 = 1.0 / g11 + SQ(ShearFactor_m * Rxy_m);
const auto g_22 = SQ(Bxy_m * hthe_m / Bpxy_m);
const auto g_33 = Rxy_m * Rxy_m;
const auto g_12 = sign_of_bp * Btxy_m * hthe_m * ShearFactor_m * Rxy_m / Bpxy_m;
const auto g_13 = ShearFactor_m * Rxy_m * Rxy_m;
const auto g_23 = sign_of_bp * Btxy_m * hthe_m * Rxy_m / Bpxy_m;

coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));

coord->setJ(hthe_m / Bpxy_m);
coord->setBxy(Bxy_m);
coord->setDx(dx_m);

return coord;
}
BoutReal ShearFactor = 1.0);

const Field2D& Rxy() const { return Rxy_m; }
const Field2D& Bpxy() const { return Bpxy_m; }
Expand Down
69 changes: 69 additions & 0 deletions src/mesh/tokamak_coordinates.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
#include <bout/tokamak_coordinates.hxx>


void TokamakCoordinates::normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) {

Rxy_m /= Lbar;
Bpxy_m /= Bbar;
Btxy_m /= Bbar;
Bxy_m /= Bbar;
hthe_m /= Lbar;
ShearFactor_m *= Lbar * Lbar * Bbar * ShearFactor;
dx_m /= Lbar * Lbar * Bbar;
}

TokamakCoordinates::TokamakCoordinates(Mesh& mesh) : mesh_m(mesh) {

mesh.get(Rxy_m, "Rxy");
// mesh->get(Zxy, "Zxy");
mesh.get(Bpxy_m, "Bpxy");
mesh.get(Btxy_m, "Btxy");
mesh.get(Bxy_m, "Bxy");
mesh.get(hthe_m, "hthe");
mesh.get(ShearFactor_m, "sinty");
mesh.get(dx_m, "dpsi");
}

BoutReal TokamakCoordinates::get_sign_of_bp() {
if (min(Bpxy_m, true) < 0.0) {
return -1.0;
}
return 1.0;
}

Coordinates* TokamakCoordinates::make_coordinates(const bool noshear, BoutReal Lbar,
BoutReal Bbar, BoutReal ShearFactor) {

normalise(Lbar, Bbar, ShearFactor);

const BoutReal sign_of_bp = get_sign_of_bp();

if (noshear) {
ShearFactor_m = 0.0;
}

auto* coord = mesh_m.getCoordinates();

const auto g11 = SQ(Rxy_m * Bpxy_m);
const auto g22 = 1.0 / SQ(hthe_m);
const auto g33 = SQ(ShearFactor_m) * g11 + SQ(Bxy_m) / g11;
const auto g12 = 0.0;
const auto g13 = -ShearFactor_m * g11;
const auto g23 = -sign_of_bp * Btxy_m / (hthe_m * Bpxy_m * Rxy_m);

const auto g_11 = 1.0 / g11 + SQ(ShearFactor_m * Rxy_m);
const auto g_22 = SQ(Bxy_m * hthe_m / Bpxy_m);
const auto g_33 = Rxy_m * Rxy_m;
const auto g_12 = sign_of_bp * Btxy_m * hthe_m * ShearFactor_m * Rxy_m / Bpxy_m;
const auto g_13 = ShearFactor_m * Rxy_m * Rxy_m;
const auto g_23 = sign_of_bp * Btxy_m * hthe_m * Rxy_m / Bpxy_m;

coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));

coord->setJ(hthe_m / Bpxy_m);
coord->setBxy(Bxy_m);
coord->setDx(dx_m);

return coord;
}

0 comments on commit 01e5816

Please sign in to comment.