From 3c4e0c8ec269b0f2dae378ff2bf8acb1bf843a9d Mon Sep 17 00:00:00 2001
From: Tom Body
Date: Wed, 19 Jun 2024 18:30:44 +0000
Subject: [PATCH 1/6] For this reason this works...
---
CMakeLists.txt | 1 +
hermes-3.cxx | 1 +
include/fieldline_geometry.hxx | 78 ++++++++++++++++++++++++++++++++++
3 files changed, 80 insertions(+)
create mode 100644 include/fieldline_geometry.hxx
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 98a96282..e2c877e4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -158,6 +158,7 @@ set(HERMES_SOURCES
include/transform.hxx
include/fixed_fraction_radiation.hxx
include/simple_pump.hxx
+ include/fieldline_geometry.hxx
)
# Compile Hermes sources into a library, which can be used for unit tests
diff --git a/hermes-3.cxx b/hermes-3.cxx
index ac5b711f..03cf12bf 100644
--- a/hermes-3.cxx
+++ b/hermes-3.cxx
@@ -77,6 +77,7 @@
#include "include/vorticity.hxx"
#include "include/zero_current.hxx"
#include "include/simple_pump.hxx"
+#include "include/fieldline_geometry.hxx"
#include
#include
#include
diff --git a/include/fieldline_geometry.hxx b/include/fieldline_geometry.hxx
new file mode 100644
index 00000000..20d8caa6
--- /dev/null
+++ b/include/fieldline_geometry.hxx
@@ -0,0 +1,78 @@
+#pragma once
+#ifndef fieldline_geometry_H
+#define fieldline_geometry_H
+
+#include "component.hxx"
+#include
+
+struct FieldlineGeometry : public Component {
+
+ FieldlineGeometry(std::string name, Options& alloptions, Solver*) : name(name) {
+
+ Options& options = alloptions[name];
+
+ const auto& units = alloptions["units"];
+ Nnorm = get(units["inv_meters_cubed"]);
+ Omega_ci = 1.0 / get(units["seconds"]);
+
+ residence_time = (options["residence_time"]
+ .doc("Pumping time-constant. Units [s]")
+ .as()
+ ) * Omega_ci;
+
+ sink_shape = (options["sink_shape"]
+ .doc("Shape of pumping sink.")
+ .withDefault(Field3D(0.0))
+ );
+
+ diagnose = options["diagnose"]
+ .doc("Output additional diagnostics?")
+ .withDefault(false);
+
+ };
+
+ void transform(Options& state) override {
+
+ Field3D species_density = getNoBoundary(state["species"][name]["density"]);
+
+ pumping_sink = (sink_shape * species_density) * (-1.0 / residence_time);
+
+ add(state["species"][name]["density_source"], pumping_sink);
+
+ };
+
+ void outputVars(Options& state) override {
+ AUTO_TRACE();
+ if (diagnose) {
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_src_shape_" + name)], sink_shape,
+ {{"long_name", "simple pump source shape"},
+ {"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_sink_" + name)], pumping_sink,
+ {{"time_dimension", "t"},
+ {"units", "m^-3 / s"},
+ {"conversion", Nnorm * Omega_ci},
+ {"long_name", "simple pump source shape"},
+ {"source", "fieldline_geometry"}});
+
+ }}
+
+ private:
+ std::string name; ///< The species name
+ Field3D sink_shape;
+ Field3D pumping_sink;
+ BoutReal Nnorm;
+ BoutReal Omega_ci;
+ BoutReal residence_time;
+ bool diagnose;
+
+};
+
+namespace {
+ RegisterComponent register_fieldline_geometry("fieldline_geometry");
+}
+
+#endif // fieldline_geometry_H
From fd541a37aea7d8ac2409a798cf7f160d7157b1b5 Mon Sep 17 00:00:00 2001
From: Tom Body
Date: Thu, 20 Jun 2024 08:26:02 +0000
Subject: [PATCH 2/6] Semi useful
---
include/fieldline_geometry.hxx | 64 +++++++++++++++-------------------
1 file changed, 28 insertions(+), 36 deletions(-)
diff --git a/include/fieldline_geometry.hxx b/include/fieldline_geometry.hxx
index 20d8caa6..f1b94e6d 100644
--- a/include/fieldline_geometry.hxx
+++ b/include/fieldline_geometry.hxx
@@ -7,25 +7,19 @@
struct FieldlineGeometry : public Component {
- FieldlineGeometry(std::string name, Options& alloptions, Solver*) : name(name) {
+ FieldlineGeometry(std::string, Options& options, Solver*) {
- Options& options = alloptions[name];
+ Options& geo_options = options["fieldline_geometry"];
- const auto& units = alloptions["units"];
- Nnorm = get(units["inv_meters_cubed"]);
- Omega_ci = 1.0 / get(units["seconds"]);
+ const auto& units = options["units"];
+ Lnorm = get(units["meters"]);
- residence_time = (options["residence_time"]
- .doc("Pumping time-constant. Units [s]")
- .as()
- ) * Omega_ci;
+ parallel_length = (geo_options["parallel_length"]
+ .doc("Parallel length as a function of grid index [m]")
+ .withDefault(Field3D(0.0))
+ ) / Lnorm;
- sink_shape = (options["sink_shape"]
- .doc("Shape of pumping sink.")
- .withDefault(Field3D(0.0))
- );
-
- diagnose = options["diagnose"]
+ diagnose = geo_options["diagnose"]
.doc("Output additional diagnostics?")
.withDefault(false);
@@ -33,12 +27,6 @@ struct FieldlineGeometry : public Component {
void transform(Options& state) override {
- Field3D species_density = getNoBoundary(state["species"][name]["density"]);
-
- pumping_sink = (sink_shape * species_density) * (-1.0 / residence_time);
-
- add(state["species"][name]["density_source"], pumping_sink);
-
};
void outputVars(Options& state) override {
@@ -46,29 +34,33 @@ struct FieldlineGeometry : public Component {
if (diagnose) {
set_with_attrs(
- state[std::string("fieldline_geometry_src_shape_" + name)], sink_shape,
- {{"long_name", "simple pump source shape"},
+ state[std::string("fieldline_geometry_parallel_length")], parallel_length,
+ {{"units", "m"},
+ {"conversion", Lnorm},
+ {"long_name", "Parallel length"},
{"source", "fieldline_geometry"}});
+
+ // set_with_attrs(
+ // state[std::string("fieldline_geometry_src_shape_" + name)], sink_shape,
+ // {{"long_name", "simple pump source shape"},
+ // {"source", "fieldline_geometry"}});
- set_with_attrs(
- state[std::string("fieldline_geometry_sink_" + name)], pumping_sink,
- {{"time_dimension", "t"},
- {"units", "m^-3 / s"},
- {"conversion", Nnorm * Omega_ci},
- {"long_name", "simple pump source shape"},
- {"source", "fieldline_geometry"}});
+ // set_with_attrs(
+ // state[std::string("fieldline_geometry_sink_" + name)], pumping_sink,
+ // {{"time_dimension", "t"},
+ // {"units", "m^-3 / s"},
+ // {"conversion", Nnorm * Omega_ci},
+ // {"long_name", "simple pump source shape"},
+ // {"source", "fieldline_geometry"}});
}}
private:
- std::string name; ///< The species name
- Field3D sink_shape;
- Field3D pumping_sink;
- BoutReal Nnorm;
- BoutReal Omega_ci;
- BoutReal residence_time;
+ BoutReal Lnorm;
bool diagnose;
+ Field3D parallel_length;
+
};
namespace {
From 5e099aefed9875aef7f2b337ec7090f646318bc7 Mon Sep 17 00:00:00 2001
From: Tom Body
Date: Thu, 20 Jun 2024 13:38:43 +0000
Subject: [PATCH 3/6] Playing around with generator for lpar
---
CMakeLists.txt | 1 +
include/fieldline_geometry.hxx | 53 +++-----------------------
src/fieldline_geometry.cxx | 68 ++++++++++++++++++++++++++++++++++
3 files changed, 75 insertions(+), 47 deletions(-)
create mode 100644 src/fieldline_geometry.cxx
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e2c877e4..6465d065 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -95,6 +95,7 @@ set(HERMES_SOURCES
src/upstream_density_feedback.cxx
src/temperature_feedback.cxx
src/detachment_controller.cxx
+ src/fieldline_geometry.cxx
src/thermal_force.cxx
src/ion_viscosity.cxx
src/transform.cxx
diff --git a/include/fieldline_geometry.hxx b/include/fieldline_geometry.hxx
index f1b94e6d..2fc355b9 100644
--- a/include/fieldline_geometry.hxx
+++ b/include/fieldline_geometry.hxx
@@ -7,59 +7,18 @@
struct FieldlineGeometry : public Component {
- FieldlineGeometry(std::string, Options& options, Solver*) {
+ FieldlineGeometry(std::string, Options& options, Solver*);
- Options& geo_options = options["fieldline_geometry"];
-
- const auto& units = options["units"];
- Lnorm = get(units["meters"]);
-
- parallel_length = (geo_options["parallel_length"]
- .doc("Parallel length as a function of grid index [m]")
- .withDefault(Field3D(0.0))
- ) / Lnorm;
-
- diagnose = geo_options["diagnose"]
- .doc("Output additional diagnostics?")
- .withDefault(false);
-
- };
-
- void transform(Options& state) override {
-
- };
-
- void outputVars(Options& state) override {
- AUTO_TRACE();
- if (diagnose) {
-
- set_with_attrs(
- state[std::string("fieldline_geometry_parallel_length")], parallel_length,
- {{"units", "m"},
- {"conversion", Lnorm},
- {"long_name", "Parallel length"},
- {"source", "fieldline_geometry"}});
-
- // set_with_attrs(
- // state[std::string("fieldline_geometry_src_shape_" + name)], sink_shape,
- // {{"long_name", "simple pump source shape"},
- // {"source", "fieldline_geometry"}});
-
- // set_with_attrs(
- // state[std::string("fieldline_geometry_sink_" + name)], pumping_sink,
- // {{"time_dimension", "t"},
- // {"units", "m^-3 / s"},
- // {"conversion", Nnorm * Omega_ci},
- // {"long_name", "simple pump source shape"},
- // {"source", "fieldline_geometry"}});
-
- }}
+ void transform(Options& state) override;
+ void outputVars(Options& state) override;
private:
BoutReal Lnorm;
bool diagnose;
- Field3D parallel_length;
+ Field3D lpar;
+ FieldGeneratorPtr B_poloidal_function;
+ Field3D B_poloidal;
};
diff --git a/src/fieldline_geometry.cxx b/src/fieldline_geometry.cxx
new file mode 100644
index 00000000..e0ba1eb7
--- /dev/null
+++ b/src/fieldline_geometry.cxx
@@ -0,0 +1,68 @@
+#include "../include/fieldline_geometry.hxx"
+#include
+#include
+#include
+#include
+#include // For outputting to log file
+
+using bout::globals::mesh;
+
+FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
+
+ Options& geo_options = options["fieldline_geometry"];
+
+ const auto& units = options["units"];
+ Lnorm = get(units["meters"]);
+
+ // lpar = (geo_options["parallel_length"]
+ // .doc("Parallel length as a function of grid index [m]")
+ // .withDefault(Field3D(0.0))
+ // ) / Lnorm;
+
+ const int MYPE = BoutComm::rank(); // Current rank
+ const int NPES = BoutComm::size(); // Number of procs
+ const int NYPE = NPES / mesh->NXPE; // Number of procs in Y
+ Coordinates *coord = mesh->getCoordinates();
+ lpar = 0;
+ auto dy = coord->dy;
+ lpar(0,0,0) = 0.5 * dy(0,0,0);
+ for (int id = 0; id <= NYPE-1; ++id) { // Iterate through each proc
+ for (int j = 0; j <= mesh->LocalNy; ++j) {
+ if (j > 0) {
+ lpar(0, j, 0) = lpar(0, j-1, 0) + 0.5 * dy(0, j-1, 0) + 0.5 * dy(0, j, 0);
+ }
+ }
+ mesh->communicate(lpar); // When proc is done, communicate to other ones
+ }
+
+ auto str = geo_options["B_poloidal"]
+ .doc("Function for B_poloidal(lpar) [T].")
+ .as();
+ B_poloidal_function = FieldFactory::get()->parse(str, &geo_options);
+
+ for (int j = mesh->ystart; j <= mesh->yend; ++j) {
+ B_poloidal(0, j, 0) = B_poloidal_function ->generate(bout::generator::Context().set("lpar", lpar(0, j, 0)));
+ }
+
+ diagnose = geo_options["diagnose"]
+ .doc("Output additional diagnostics?")
+ .withDefault(false);
+
+}
+
+void FieldlineGeometry::transform(Options& state) {
+
+}
+
+void FieldlineGeometry::outputVars(Options& state) {
+ AUTO_TRACE();
+ if (diagnose) {
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_parallel_length")], lpar,
+ {{"units", "m"},
+ // {"conversion", Lnorm},
+ {"long_name", "Parallel length"},
+ {"source", "fieldline_geometry"}});
+
+}}
\ No newline at end of file
From 58840c4524fedf2c9e80d7f9728085445f510367 Mon Sep 17 00:00:00 2001
From: Tom Body
Date: Thu, 20 Jun 2024 14:26:42 +0000
Subject: [PATCH 4/6] Correctly parse lpar from input file
---
include/fieldline_geometry.hxx | 4 ++--
src/fieldline_geometry.cxx | 20 +++++++++++---------
2 files changed, 13 insertions(+), 11 deletions(-)
diff --git a/include/fieldline_geometry.hxx b/include/fieldline_geometry.hxx
index 2fc355b9..1d1134b5 100644
--- a/include/fieldline_geometry.hxx
+++ b/include/fieldline_geometry.hxx
@@ -16,9 +16,9 @@ struct FieldlineGeometry : public Component {
BoutReal Lnorm;
bool diagnose;
- Field3D lpar;
+ Field3D lpar{0.0};
FieldGeneratorPtr B_poloidal_function;
- Field3D B_poloidal;
+ Field3D B_poloidal{0.0};
};
diff --git a/src/fieldline_geometry.cxx b/src/fieldline_geometry.cxx
index e0ba1eb7..62df3a77 100644
--- a/src/fieldline_geometry.cxx
+++ b/src/fieldline_geometry.cxx
@@ -14,11 +14,6 @@ FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
const auto& units = options["units"];
Lnorm = get(units["meters"]);
- // lpar = (geo_options["parallel_length"]
- // .doc("Parallel length as a function of grid index [m]")
- // .withDefault(Field3D(0.0))
- // ) / Lnorm;
-
const int MYPE = BoutComm::rank(); // Current rank
const int NPES = BoutComm::size(); // Number of procs
const int NYPE = NPES / mesh->NXPE; // Number of procs in Y
@@ -36,12 +31,14 @@ FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
}
auto str = geo_options["B_poloidal"]
- .doc("Function for B_poloidal(lpar) [T].")
+ .doc("Function for B_poloidal(lpar) (for B_poloidal in [T] and lpar in [m]).")
.as();
B_poloidal_function = FieldFactory::get()->parse(str, &geo_options);
- for (int j = mesh->ystart; j <= mesh->yend; ++j) {
- B_poloidal(0, j, 0) = B_poloidal_function ->generate(bout::generator::Context().set("lpar", lpar(0, j, 0)));
+ B_poloidal.allocate();
+
+ BOUT_FOR(i, B_poloidal.getRegion("RGN_ALL")) {
+ B_poloidal[i] = B_poloidal_function->generate(bout::generator::Context().set("lpar", lpar[i]));
}
diagnose = geo_options["diagnose"]
@@ -61,8 +58,13 @@ void FieldlineGeometry::outputVars(Options& state) {
set_with_attrs(
state[std::string("fieldline_geometry_parallel_length")], lpar,
{{"units", "m"},
- // {"conversion", Lnorm},
{"long_name", "Parallel length"},
{"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_B_poloidal")], B_poloidal,
+ {{"units", "T"},
+ {"long_name", "B poloidal"},
+ {"source", "fieldline_geometry"}});
}}
\ No newline at end of file
From cf10bedec839b3d72616371431496fa1cd340e47 Mon Sep 17 00:00:00 2001
From: Tom Body
Date: Thu, 20 Jun 2024 14:30:38 +0000
Subject: [PATCH 5/6] Fix computation of lpar
---
src/fieldline_geometry.cxx | 9 ++++++++-
1 file changed, 8 insertions(+), 1 deletion(-)
diff --git a/src/fieldline_geometry.cxx b/src/fieldline_geometry.cxx
index 62df3a77..83ecc1be 100644
--- a/src/fieldline_geometry.cxx
+++ b/src/fieldline_geometry.cxx
@@ -14,11 +14,13 @@ FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
const auto& units = options["units"];
Lnorm = get(units["meters"]);
+ // Get lpar
const int MYPE = BoutComm::rank(); // Current rank
const int NPES = BoutComm::size(); // Number of procs
const int NYPE = NPES / mesh->NXPE; // Number of procs in Y
Coordinates *coord = mesh->getCoordinates();
lpar = 0;
+ BoutReal offset = 0; // Offset to ensure ylow domain boundary starts at 0
auto dy = coord->dy;
lpar(0,0,0) = 0.5 * dy(0,0,0);
for (int id = 0; id <= NYPE-1; ++id) { // Iterate through each proc
@@ -27,8 +29,13 @@ FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
lpar(0, j, 0) = lpar(0, j-1, 0) + 0.5 * dy(0, j-1, 0) + 0.5 * dy(0, j, 0);
}
}
- mesh->communicate(lpar); // When proc is done, communicate to other ones
+ mesh->communicate(lpar); // Communicate across guard cells so other procs keep counting
+ if (MYPE == 0) {
+ offset = (lpar(0,1,0) + lpar(0,2,0)) / 2; // Offset lives on first proc
+ }
}
+ MPI_Bcast(&offset, 1, MPI_DOUBLE, 0, BoutComm::get()); // Ensure all procs get offset
+ lpar -= offset;
auto str = geo_options["B_poloidal"]
.doc("Function for B_poloidal(lpar) (for B_poloidal in [T] and lpar in [m]).")
From 80084d0004c461ba0a31711e13487e647aebe794 Mon Sep 17 00:00:00 2001
From: Tom Body
Date: Thu, 20 Jun 2024 17:27:06 +0000
Subject: [PATCH 6/6] Calculate geometry factors according to Derks 2024
---
include/fieldline_geometry.hxx | 11 +++-
src/fieldline_geometry.cxx | 104 +++++++++++++++++++++++++++++----
2 files changed, 103 insertions(+), 12 deletions(-)
diff --git a/include/fieldline_geometry.hxx b/include/fieldline_geometry.hxx
index 1d1134b5..270d33d8 100644
--- a/include/fieldline_geometry.hxx
+++ b/include/fieldline_geometry.hxx
@@ -17,8 +17,17 @@ struct FieldlineGeometry : public Component {
bool diagnose;
Field3D lpar{0.0};
- FieldGeneratorPtr B_poloidal_function;
+
Field3D B_poloidal{0.0};
+ Field3D B_toroidal{0.0};
+ Field3D major_radius{0.0};
+ Field3D effective_flux_exp{0.0};
+ Field3D poloidal_SOL_width{0.0};
+ BoutReal lambda_q_omp;
+
+ Field3D B_total{0.0};
+ Field3D pitch_angle{0.0};
+ Field3D area_external{0.0};
};
diff --git a/src/fieldline_geometry.cxx b/src/fieldline_geometry.cxx
index 83ecc1be..95d340d2 100644
--- a/src/fieldline_geometry.cxx
+++ b/src/fieldline_geometry.cxx
@@ -7,19 +7,14 @@
using bout::globals::mesh;
-FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
-
- Options& geo_options = options["fieldline_geometry"];
-
- const auto& units = options["units"];
- Lnorm = get(units["meters"]);
-
+Field3D calculateLpar() {
// Get lpar
const int MYPE = BoutComm::rank(); // Current rank
const int NPES = BoutComm::size(); // Number of procs
const int NYPE = NPES / mesh->NXPE; // Number of procs in Y
Coordinates *coord = mesh->getCoordinates();
- lpar = 0;
+ Field3D lpar{0.0};
+
BoutReal offset = 0; // Offset to ensure ylow domain boundary starts at 0
auto dy = coord->dy;
lpar(0,0,0) = 0.5 * dy(0,0,0);
@@ -37,17 +32,68 @@ FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
MPI_Bcast(&offset, 1, MPI_DOUBLE, 0, BoutComm::get()); // Ensure all procs get offset
lpar -= offset;
- auto str = geo_options["B_poloidal"]
+ return lpar;
+}
+
+FieldlineGeometry::FieldlineGeometry(std::string, Options& options, Solver*) {
+
+ Options& geo_options = options["fieldline_geometry"];
+ Coordinates *coord = mesh->getCoordinates();
+
+ const auto& units = options["units"];
+ Lnorm = get(units["meters"]);
+
+ lpar = calculateLpar();
+
+ lambda_q_omp = geo_options["lambda_q_omp"]
+ .doc("Heat-flux decay length at outboard midplane [m]")
+ .as();
+
+ auto B_poloidal_str = geo_options["B_poloidal"]
.doc("Function for B_poloidal(lpar) (for B_poloidal in [T] and lpar in [m]).")
.as();
- B_poloidal_function = FieldFactory::get()->parse(str, &geo_options);
+ FieldGeneratorPtr B_poloidal_function = FieldFactory::get()->parse(B_poloidal_str, &geo_options);
+
+ auto B_toroidal_str = geo_options["B_toroidal"]
+ .doc("Function for B_toroidal(lpar) (for B_toroidal in [T] and lpar in [m]).")
+ .as();
+ FieldGeneratorPtr B_toroidal_function = FieldFactory::get()->parse(B_toroidal_str, &geo_options);
+
+ auto major_radius_str = geo_options["major_radius"]
+ .doc("Function for major_radius(lpar) (for major_radius in [m] and lpar in [m]).")
+ .as();
+ FieldGeneratorPtr major_radius_function = FieldFactory::get()->parse(major_radius_str, &geo_options);
+
+ auto effective_flux_exp_str = geo_options["effective_flux_exp"]
+ .doc("Function for effective_flux_exp(lpar) (for effective_flux_exp in [~] and lpar in [m]).")
+ .as();
+ FieldGeneratorPtr effective_flux_exp_function = FieldFactory::get()->parse(effective_flux_exp_str, &geo_options);
B_poloidal.allocate();
+ B_toroidal.allocate();
+ major_radius.allocate();
+ effective_flux_exp.allocate();
- BOUT_FOR(i, B_poloidal.getRegion("RGN_ALL")) {
+ BOUT_FOR(i, lpar.getRegion("RGN_ALL")) {
B_poloidal[i] = B_poloidal_function->generate(bout::generator::Context().set("lpar", lpar[i]));
+ B_toroidal[i] = B_toroidal_function->generate(bout::generator::Context().set("lpar", lpar[i]));
+ major_radius[i] = major_radius_function->generate(bout::generator::Context().set("lpar", lpar[i]));
+ effective_flux_exp[i] = effective_flux_exp_function->generate(bout::generator::Context().set("lpar", lpar[i]));
}
+ B_total = sqrt(B_toroidal*B_toroidal + B_poloidal*B_poloidal);
+ pitch_angle = B_poloidal / B_total;
+
+ auto dy = coord->dy;
+ area_external = pitch_angle * dy * 2.0 * PI * major_radius;
+
+ BoutReal upstream_pitch_angle = pitch_angle(0, mesh->ystart, 0);
+ BoutReal upstream_effective_flux_exp = effective_flux_exp(0, mesh->ystart, 0);
+
+ // "lambda_q" of equation 6 of Derks 2024
+ // N.b. B_trans is the inverse of effective_flux_exp.
+ poloidal_SOL_width = lambda_q_omp * (upstream_pitch_angle / pitch_angle) * (effective_flux_exp / upstream_effective_flux_exp);
+
diagnose = geo_options["diagnose"]
.doc("Output additional diagnostics?")
.withDefault(false);
@@ -73,5 +119,41 @@ void FieldlineGeometry::outputVars(Options& state) {
{{"units", "T"},
{"long_name", "B poloidal"},
{"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_B_toroidal")], B_toroidal,
+ {{"units", "T"},
+ {"long_name", "B toroidal"},
+ {"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_major_radius")], major_radius,
+ {{"units", "m"},
+ {"long_name", "major radius"},
+ {"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_effective_flux_exp")], effective_flux_exp,
+ {{"units", ""},
+ {"long_name", "effective flux expansion due to cross-field transport"},
+ {"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_pitch_angle")], pitch_angle,
+ {{"units", ""},
+ {"long_name", "B-poloidal / B-total"},
+ {"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_area_external")], area_external,
+ {{"units", "m^2"},
+ {"long_name", "area perpendicular to poloidal SOL width"},
+ {"source", "fieldline_geometry"}});
+
+ set_with_attrs(
+ state[std::string("fieldline_geometry_poloidal_SOL_width")], poloidal_SOL_width,
+ {{"units", "m"},
+ {"long_name", "poloidal scrape-off-layer width"},
+ {"source", "fieldline_geometry"}});
}}
\ No newline at end of file