Skip to content

Commit

Permalink
Merge branch 'always-normalise' into dry-up-examples
Browse files Browse the repository at this point in the history
  • Loading branch information
@tomchapman
tomchapman committed Dec 16, 2024
2 parents ae7f55f + 8292377 commit 72ce52d
Show file tree
Hide file tree
Showing 15 changed files with 178 additions and 178 deletions.
99 changes: 51 additions & 48 deletions examples/6field-simple/elm_6f.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,7 +366,7 @@ class Elm_6f : public PhysicsModel {
result = Grad_par(f, loc);

if (nonlinear) {
result -= bracket(Psi, f, bm_mag) * tokamak_coordinates_factory.get_Bxy();
result -= bracket(Psi, f, bm_mag) * tokamak_coordinates_factory.Bxy();
}
}

Expand Down Expand Up @@ -678,7 +678,7 @@ class Elm_6f : public PhysicsModel {

if (noshear) {
if (include_curvature) {
b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
}
}

Expand All @@ -688,7 +688,7 @@ class Elm_6f : public PhysicsModel {
if (not mesh->IncIntShear) {
// Dimits style, using local coordinate system
if (include_curvature) {
b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
}
}

Expand Down Expand Up @@ -832,11 +832,9 @@ class Elm_6f : public PhysicsModel {
dump.add(sp_length, "sp_length", 1);
}

J0 = SI::mu0 * Lbar * J0 / tokamak_coordinates_factory.get_Bxy();
J0 = SI::mu0 * Lbar * J0 / tokamak_coordinates_factory.Bxy();
P0 = P0 / (SI::kb * (Tibar + Tebar) * eV_K / 2. * Nbar * density);

tokamak_coordinates_factory.normalise(Lbar, Bbar);

b0xcv.x /= Bbar;
b0xcv.y *= Lbar * Lbar;
b0xcv.z *= Lbar * Lbar;
Expand Down Expand Up @@ -906,7 +904,10 @@ class Elm_6f : public PhysicsModel {
q95 = q95_input; // use a constant for test
} else {
if (local_q) {
q95 = abs(tokamak_coordinates_factory.get_hthe() * tokamak_coordinates_factory.get_Btxy() / tokamak_coordinates_factory.get_Bpxy()) * q_alpha;
q95 = abs(tokamak_coordinates_factory.hthe()
* tokamak_coordinates_factory.Btxy()
/ tokamak_coordinates_factory.Bpxy())
* q_alpha;
} else {
output.write("\tUsing q profile from grid.\n");
if (mesh->get(q95, "q")) {
Expand Down Expand Up @@ -1028,36 +1029,37 @@ class Elm_6f : public PhysicsModel {
}

/**************** CALCULATE METRICS ******************/
const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);

const auto& coord =
tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lbar, Bbar);

//////////////////////////////////////////////////////////////
// SHIFTED RADIAL COORDINATES

if (mesh->IncIntShear) {
// BOUT-06 style, using d/dx = d/dpsi + I * d/dz
coord->setIntShiftTorsion(tokamak_coordinates_factory.get_ShearFactor());
coord->setIntShiftTorsion(tokamak_coordinates_factory.ShearFactor());
}

// Set B field vector

B0vec.covariant = false;
B0vec.x = 0.;
B0vec.y = tokamak_coordinates_factory.get_Bpxy() / tokamak_coordinates_factory.get_hthe();
B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe();
B0vec.z = 0.;

// Set V0vec field vector

V0vec.covariant = false;
V0vec.x = 0.;
V0vec.y = Vp0 / tokamak_coordinates_factory.get_hthe();
V0vec.z = Vt0 / tokamak_coordinates_factory.get_Rxy();
V0vec.y = Vp0 / tokamak_coordinates_factory.hthe();
V0vec.z = Vt0 / tokamak_coordinates_factory.Rxy();

// Set V0eff field vector

V0eff.covariant = false;
V0eff.x = 0.;
V0eff.y = -(tokamak_coordinates_factory.get_Btxy() / (tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy())) * (Vp0 * tokamak_coordinates_factory.get_Btxy() - Vt0 * tokamak_coordinates_factory.get_Bpxy()) / tokamak_coordinates_factory.get_hthe();
V0eff.z = (tokamak_coordinates_factory.get_Bpxy() / (tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy())) * (Vp0 * tokamak_coordinates_factory.get_Btxy() - Vt0 * tokamak_coordinates_factory.get_Bpxy()) / tokamak_coordinates_factory.get_Rxy();
V0eff.y = -(tokamak_coordinates_factory.Btxy() / (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy())) * (Vp0 * tokamak_coordinates_factory.Btxy() - Vt0 * tokamak_coordinates_factory.Bpxy()) / tokamak_coordinates_factory.hthe();
V0eff.z = (tokamak_coordinates_factory.Bpxy() / (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy())) * (Vp0 * tokamak_coordinates_factory.Btxy() - Vt0 * tokamak_coordinates_factory.Bpxy()) / tokamak_coordinates_factory.Rxy();

Pe.setBoundary("P");
Pi.setBoundary("P");
Expand Down Expand Up @@ -1099,10 +1101,10 @@ class Elm_6f : public PhysicsModel {
if (diamag && diamag_phi0) {
if (experiment_Er) { // get phi0 from grid file
mesh->get(phi0, "Phi_0");
phi0 /= tokamak_coordinates_factory.get_Bxy() * Lbar * Va;
phi0 /= tokamak_coordinates_factory.Bxy() * Lbar * Va;
} else {
// Stationary equilibrium plasma. ExB velocity balances diamagnetic drift
phi0 = -Upara0 * Pi0 / tokamak_coordinates_factory.get_Bxy() / N0;
phi0 = -Upara0 * Pi0 / tokamak_coordinates_factory.Bxy() / N0;
}
SAVE_ONCE(phi0);
}
Expand All @@ -1112,7 +1114,7 @@ class Elm_6f : public PhysicsModel {
SAVE_ONCE(J0, P0);
SAVE_ONCE(density, Lbar, Bbar, Tbar);
SAVE_ONCE(Tibar, Tebar, Nbar);
Field2D tmp = tokamak_coordinates_factory.get_Bxy();
Field2D tmp = tokamak_coordinates_factory.Bxy();
SAVE_ONCE(Va, tmp);
SAVE_ONCE(Ti0, Te0, N0);

Expand All @@ -1136,13 +1138,13 @@ class Elm_6f : public PhysicsModel {
Field2D logn0 = laplace_alpha * N0;
Field3D Ntemp;
Ntemp = N0;
ubyn = U * tokamak_coordinates_factory.get_Bxy() / Ntemp;
ubyn = U * tokamak_coordinates_factory.Bxy() / Ntemp;
// Phi should be consistent with U
if (laplace_alpha <= 0.0) {
phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.get_Bxy();
phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy();
} else {
phiSolver->setCoefC(logn0);
phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.get_Bxy();
phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy();
}
}

Expand Down Expand Up @@ -1191,6 +1193,7 @@ class Elm_6f : public PhysicsModel {

return 0;
}

int rhs(BoutReal UNUSED(t)) override {

Coordinates* coord = mesh->getCoordinates();
Expand All @@ -1216,17 +1219,17 @@ class Elm_6f : public PhysicsModel {

// Field2D lap_temp=0.0;
Field2D logn0 = laplace_alpha * N0;
ubyn = U * tokamak_coordinates_factory.get_Bxy() / N0;
ubyn = U * tokamak_coordinates_factory.Bxy() / N0;
if (diamag) {
ubyn -= Upara0 / N0 * Delp2(Pi) / tokamak_coordinates_factory.get_Bxy();
ubyn -= Upara0 / N0 * Delp2(Pi) / tokamak_coordinates_factory.Bxy();
mesh->communicate(ubyn);
ubyn.applyBoundary();
}
// Invert laplacian for phi
if (laplace_alpha > 0.0) {
phiSolver->setCoefC(logn0);
}
phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.get_Bxy();
phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy();

mesh->communicate(phi);

Expand Down Expand Up @@ -1322,9 +1325,9 @@ class Elm_6f : public PhysicsModel {

if (compress0) {
if (nonlinear) {
Vepar = Vipar - tokamak_coordinates_factory.get_Bxy() * (Jpar) / N_tmp * Vepara;
Vepar = Vipar - tokamak_coordinates_factory.Bxy() * (Jpar) / N_tmp * Vepara;
} else {
Vepar = Vipar - tokamak_coordinates_factory.get_Bxy() * (Jpar) / N0 * Vepara;
Vepar = Vipar - tokamak_coordinates_factory.Bxy() * (Jpar) / N0 * Vepara;
Vepar.applyBoundary();
mesh->communicate(Vepar);
}
Expand Down Expand Up @@ -1362,13 +1365,13 @@ class Elm_6f : public PhysicsModel {
ddt(Psi) = 0.0;

if (spitzer_resist) {
ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.get_Bxy() * phi) / tokamak_coordinates_factory.get_Bxy() - eta_spitzer * Jpar;
ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.Bxy() * phi) / tokamak_coordinates_factory.Bxy() - eta_spitzer * Jpar;
} else {
ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.get_Bxy() * phi) / tokamak_coordinates_factory.get_Bxy() - eta * Jpar;
ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.Bxy() * phi) / tokamak_coordinates_factory.Bxy() - eta * Jpar;
}

if (diamag) {
ddt(Psi) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Psi, bm_exb); // Equilibrium flow
ddt(Psi) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Psi, bm_exb); // Equilibrium flow
}

// Hyper-resistivity
Expand All @@ -1385,18 +1388,18 @@ class Elm_6f : public PhysicsModel {

ddt(U) = 0.0;

ddt(U) = -SQ(tokamak_coordinates_factory.get_Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates_factory.get_Bxy(); // Grad j term
ddt(U) = -SQ(tokamak_coordinates_factory.Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates_factory.Bxy(); // Grad j term

ddt(U) += 2.0 * Upara1 * b0xcv * Grad(P); // curvature term

ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar); // b dot grad j
ddt(U) += SQ(tokamak_coordinates_factory.Bxy()) * Grad_parP(Jpar); // b dot grad j

if (diamag) {
ddt(U) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, U, bm_exb); // Equilibrium flow
ddt(U) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, U, bm_exb); // Equilibrium flow
}

if (nonlinear) {
ddt(U) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, U, bm_exb); // Advection
ddt(U) -= bracket(tokamak_coordinates_factory.Bxy() * phi, U, bm_exb); // Advection
}

// parallel hyper-viscous diffusion for vector potential
Expand Down Expand Up @@ -1448,18 +1451,18 @@ class Elm_6f : public PhysicsModel {

ddt(Ni) = 0.0;

ddt(Ni) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, N0, bm_exb);
ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi, N0, bm_exb);

if (diamag) {
ddt(Ni) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Ni, bm_exb); // Equilibrium flow
ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Ni, bm_exb); // Equilibrium flow
}

if (nonlinear) {
ddt(Ni) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Ni, bm_exb); // Advection
ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Ni, bm_exb); // Advection
}

if (compress0) {
ddt(Ni) -= N0 * tokamak_coordinates_factory.get_Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.get_Bxy());
ddt(Ni) -= N0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.Bxy());
}

// 4th order Parallel diffusion terms
Expand All @@ -1478,18 +1481,18 @@ class Elm_6f : public PhysicsModel {

ddt(Ti) = 0.0;

ddt(Ti) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Ti0, bm_exb);
ddt(Ti) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Ti0, bm_exb);

if (diamag) {
ddt(Ti) -= bracket(phi0 * tokamak_coordinates_factory.get_Bxy(), Ti, bm_exb); // Equilibrium flow
ddt(Ti) -= bracket(phi0 * tokamak_coordinates_factory.Bxy(), Ti, bm_exb); // Equilibrium flow
}

if (nonlinear) {
ddt(Ti) -= bracket(phi * tokamak_coordinates_factory.get_Bxy(), Ti, bm_exb); // Advection
ddt(Ti) -= bracket(phi * tokamak_coordinates_factory.Bxy(), Ti, bm_exb); // Advection
}

if (compress0) {
ddt(Ti) -= 2.0 / 3.0 * Ti0 * tokamak_coordinates_factory.get_Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.get_Bxy());
ddt(Ti) -= 2.0 / 3.0 * Ti0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.Bxy());
}

if (diffusion_par > 0.0) {
Expand All @@ -1514,18 +1517,18 @@ class Elm_6f : public PhysicsModel {

ddt(Te) = 0.0;

ddt(Te) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Te0, bm_exb);
ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Te0, bm_exb);

if (diamag) {
ddt(Te) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Te, bm_exb); // Equilibrium flow
ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Te, bm_exb); // Equilibrium flow
}

if (nonlinear) {
ddt(Te) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Te, bm_exb); // Advection
ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Te, bm_exb); // Advection
}

if (compress0) {
ddt(Te) -= 2.0 / 3.0 * Te0 * tokamak_coordinates_factory.get_Bxy() * Grad_parP(Vepar / tokamak_coordinates_factory.get_Bxy());
ddt(Te) -= 2.0 / 3.0 * Te0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vepar / tokamak_coordinates_factory.Bxy());
}

if (diffusion_par > 0.0) {
Expand All @@ -1548,14 +1551,14 @@ class Elm_6f : public PhysicsModel {
ddt(Vipar) = 0.0;

ddt(Vipar) -= Vipara * Grad_parP(P) / N0;
ddt(Vipar) += Vipara * bracket(Psi, P0, bm_mag) * tokamak_coordinates_factory.get_Bxy() / N0;
ddt(Vipar) += Vipara * bracket(Psi, P0, bm_mag) * tokamak_coordinates_factory.Bxy() / N0;

if (diamag) {
ddt(Vipar) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Vipar, bm_exb);
ddt(Vipar) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Vipar, bm_exb);
}

if (nonlinear) {
ddt(Vipar) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Vipar, bm_exb);
ddt(Vipar) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Vipar, bm_exb);
}

// parallel hyper-viscous diffusion for vector potential
Expand Down
12 changes: 6 additions & 6 deletions examples/conducting-wall-mode/cwm.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,8 +129,8 @@ class CWM : public PhysicsModel {
Te0 /= Te_x;

// Normalise geometry
tokamak_coordinates_factory.normalise(rho_s, bmag / 1e4, ShearFactor);
coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, rho_s,
bmag / 1e4, ShearFactor);

// Set nu
nu = nu_hat * Ni0 / pow(Te0, 1.5);
Expand All @@ -141,10 +141,10 @@ class CWM : public PhysicsModel {
// add evolving variables to the communication object
SOLVE_FOR(rho, te);

Field2D Rxy = tokamak_coordinates_factory.get_Rxy();
Field2D Bpxy = tokamak_coordinates_factory.get_Bpxy();
Field2D Btxy = tokamak_coordinates_factory.get_Btxy();
Field2D hthe = tokamak_coordinates_factory.get_hthe();
Field2D Rxy = tokamak_coordinates_factory.Rxy();
Field2D Bpxy = tokamak_coordinates_factory.Bpxy();
Field2D Btxy = tokamak_coordinates_factory.Btxy();
Field2D hthe = tokamak_coordinates_factory.hthe();
SAVE_ONCE(Rxy, Bpxy, Btxy, Zxy, hthe);
SAVE_ONCE(nu_hat, hthe0);

Expand Down
4 changes: 2 additions & 2 deletions examples/constraints/alfven-wave/alfven.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,8 +171,8 @@ class Alfven : public PhysicsModel {
}

auto tokamak_coordinates_factory = TokamakCoordinatesFactory(*mesh);
tokamak_coordinates_factory.normalise(Lnorm, Bnorm);
const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
const auto& coord =
tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lnorm, Bnorm);
}
};

Expand Down
6 changes: 3 additions & 3 deletions examples/dalf3/dalf3.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ class DALF3 : public PhysicsModel {

if (lowercase(ptstr) == "shifted") {
// Dimits style, using local coordinate system
b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
noshear = true;
}

Expand Down Expand Up @@ -222,8 +222,8 @@ class DALF3 : public PhysicsModel {
b0xcv.z *= rho_s * rho_s;

// Metrics
tokamak_coordinates_factory.normalise(rho_s, Bnorm);
const auto& coord = tokamak_coordinates_factory.make_tokamak_coordinates(noshear);
const auto& coord =
tokamak_coordinates_factory.make_tokamak_coordinates(noshear, rho_s, Bnorm);

SOLVE_FOR3(Vort, Pe, Vpar);
comms.add(Vort, Pe, Vpar);
Expand Down
Loading

0 comments on commit 72ce52d

Please sign in to comment.