Skip to content

Non-Boussinesq vorticity #54

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

Closed
johnomotani wants to merge 1 commit into from
Closed

Non-Boussinesq vorticity #54

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
4 changes: 4 additions & 0 deletions include/vorticity.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ private:
bool exb_advection; // Include nonlinear ExB advection?
bool diamagnetic; // Include diamagnetic current?
bool diamagnetic_polarisation; // Include diamagnetic drift in polarisation current
bool boussinesq; // Use 'Boussinesq' approximation?
BoutReal average_atomic_mass; // Weighted average atomic mass, for polarisaion current (Boussinesq approximation)
bool poloidal_flows; ///< Include poloidal ExB flow?
bool bndry_flux; ///< Allow flows through radial boundaries?
Expand All @@ -83,6 +84,9 @@ private:
Vector2D Curlb_B; // Curvature vector Curl(b/B)
BoutReal hyper_z; ///< Hyper-viscosity in Z

// Intermediate variables to carry between transform() and finally()
Field3D N_sum, Pi_sum;

// Diagnostic outputs
Field3D DivJdia, DivJcol; // Divergence of diamagnetic and collisional current
};
Expand Down
219 changes: 158 additions & 61 deletions src/vorticity.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,11 @@ Vorticity::Vorticity(std::string name, Options &alloptions, Solver *solver) {
.doc("Include diamagnetic drift in polarisation current?")
.withDefault<bool>(true);

boussinesq =
options["boussinesq"]
.doc("Use 'Boussinesq' approximation?")
.withDefault<bool>(true);

collisional_friction =
options["collisional_friction"]
.doc("Damp vorticity based on mass-weighted collision frequency")
Expand Down Expand Up @@ -99,8 +104,15 @@ Vorticity::Vorticity(std::string name, Options &alloptions, Solver *solver) {
laplacexy->setCoefs(average_atomic_mass / SQ(coord->Bxy), 0.0);
}
phiSolver = Laplacian::create(&options["laplacian"]);
// Set coefficients for Boussinesq solve
phiSolver->setCoefC(average_atomic_mass / SQ(coord->Bxy));
if (boussinesq) {
// Set coefficients for Boussinesq solve
phiSolver->setCoefC(average_atomic_mass / SQ(coord->Bxy));
} else {
if (not phiSolver->uses3DCoefs()) {
throw BoutException("Non-Boussinesq solve requires phiSolver to support 3D "
"coefficients.");
}
}

if (phi_boundary_relax) {
// Set the last update time to -1, so it will reset
Expand Down Expand Up @@ -175,6 +187,7 @@ void Vorticity::transform(Options &state) {
AUTO_TRACE();

auto& fields = state["fields"];
auto& allspecies = state["species"];

// Set the boundary of phi. Both 2D and 3D fields are kept, though the 3D field
// is constant in Z. This is for efficiency, to reduce the number of conversions.
Expand All @@ -184,7 +197,6 @@ void Vorticity::transform(Options &state) {
if (diamagnetic_polarisation) {
// Diamagnetic term in vorticity. Note this is weighted by the mass
// This includes all species, including electrons
Options& allspecies = state["species"];
for (auto& kv : allspecies.getChildren()) {
Options& species = allspecies[kv.first]; // Note: need non-const

Expand Down Expand Up @@ -336,60 +348,121 @@ void Vorticity::transform(Options &state) {
}
}

// Update boundary conditions. Two issues:
// 1) Solving here for phi + Pi, and then subtracting Pi from the result
// The boundary values should therefore include Pi
// 2) The INVERT_SET flag takes the value in the guard (boundary) cell
// and sets the boundary between cells to this value.
// This shift by 1/2 grid cell is important.
if (boussinesq) {
// Update boundary conditions. Two issues:
// 1) Solving here for phi + Pi, and then subtracting Pi from the result
// The boundary values should therefore include Pi
// 2) The INVERT_SET flag takes the value in the guard (boundary) cell
// and sets the boundary between cells to this value.
// This shift by 1/2 grid cell is important.

Field3D phi_plus_pi = phi + Pi_sum;
Field3D phi_plus_pi = phi + Pi_sum;

if (mesh->firstX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
// Average phi + Pi at the boundary, and set the boundary cell
// to this value. The phi solver will then put the value back
// onto the cell mid-point
phi_plus_pi(mesh->xstart - 1, j, k) =
0.5
* (phi_plus_pi(mesh->xstart - 1, j, k) +
phi_plus_pi(mesh->xstart, j, k));
if (mesh->firstX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
// Average phi + Pi at the boundary, and set the boundary cell
// to this value. The phi solver will then put the value back
// onto the cell mid-point
phi_plus_pi(mesh->xstart - 1, j, k) =
0.5
* (phi_plus_pi(mesh->xstart - 1, j, k) +
phi_plus_pi(mesh->xstart, j, k));
}
}
}
}

if (mesh->lastX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
phi_plus_pi(mesh->xend + 1, j, k) =
0.5
* (phi_plus_pi(mesh->xend + 1, j, k) +
phi_plus_pi(mesh->xend, j, k));
if (mesh->lastX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
phi_plus_pi(mesh->xend + 1, j, k) =
0.5
* (phi_plus_pi(mesh->xend + 1, j, k) +
phi_plus_pi(mesh->xend, j, k));
}
}
}
}

// Calculate potential
if (split_n0) {
////////////////////////////////////////////
// Split into axisymmetric and non-axisymmetric components
Field2D Vort2D = DC(Vort); // n=0 component
Field2D phi_plus_pi_2d = DC(phi_plus_pi);
phi_plus_pi -= phi_plus_pi_2d;
// Calculate potential
if (split_n0) {
////////////////////////////////////////////
// Split into axisymmetric and non-axisymmetric components
Field2D Vort2D = DC(Vort); // n=0 component
Field2D phi_plus_pi_2d = DC(phi_plus_pi);
phi_plus_pi -= phi_plus_pi_2d;

phi_plus_pi_2d = laplacexy->solve(Vort2D, phi_plus_pi_2d);
phi_plus_pi_2d = laplacexy->solve(Vort2D, phi_plus_pi_2d);

// Solve non-axisymmetric part using X-Z solver
phi = phi_plus_pi_2d
+ phiSolver->solve((Vort - Vort2D) * (Bsq / average_atomic_mass),
phi_plus_pi)
- Pi_sum;
// Solve non-axisymmetric part using X-Z solver
phi = phi_plus_pi_2d
+ phiSolver->solve((Vort - Vort2D) * (Bsq / average_atomic_mass), phi_plus_pi)
- Pi_sum;

} else {
phi = phiSolver->solve(Vort * (Bsq / average_atomic_mass),
phi_plus_pi)
- Pi_sum;
}
} else {
phi = phiSolver->solve(Vort * (Bsq / average_atomic_mass),
phi_plus_pi)
- Pi_sum;
// If non-Boussinesq, subtract Pi contribution to Vort first, then solve for phi, not
// phi_plus_pi.
// Need to do this way (subtract Div(B^-2 Grad_perp(pi)) from rhs) for non-Boussinesq
// solve, as coefficients of Grad_perp(phi) and Grad_perp(Pi) terms are not the same
// in non-Boussinesq vorticity

if (split_n0) {
throw BoutException("split_n0 not implemented (yet?) for non-Boussinesq case");
}

// Update boundary conditions.
// The INVERT_SET flag takes the value in the guard (boundary) cell and sets the
// boundary between cells to this value. This shift by 1/2 grid cell is
// important.

Field3D phi_boundary = copy(phi);
if (mesh->firstX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
// Average phi + Pi at the boundary, and set the boundary cell
// to this value. The phi solver will then put the value back
// onto the cell mid-point
phi_boundary(mesh->xstart - 1, j, k) =
0.5
* (phi_boundary(mesh->xstart - 1, j, k) +
phi_boundary(mesh->xstart, j, k));
}
}
}

if (mesh->lastX()) {
for (int j = mesh->ystart; j <= mesh->yend; j++) {
for (int k = 0; k < mesh->LocalNz; k++) {
phi_boundary(mesh->xend + 1, j, k) =
0.5
* (phi_boundary(mesh->xend + 1, j, k) +
phi_boundary(mesh->xend, j, k));
}
}
}

N_sum = 0.0;
for (auto& kv : allspecies.getChildren()) {
Options& species = allspecies[kv.first]; // Note: need non-const

if (!(IS_SET_NOBOUNDARY(species["density"]) and
species.isSet("charge") and species.isSet("AA"))) {
continue; // No density, charge or mass -> no polarisation current
}
auto N = GET_NOBOUNDARY(Field3D, species["density"]);
N_sum += N;
}
N_sum.applyBoundary("neumann");
auto n_B2 = N_sum/Bsq;
phiSolver->setCoefC(n_B2); // Set when initialised
// phi_boundary here is equivalent to withBoundary(phi, phi_boundary) because
// of the way phi_boundary is initialised
phi = phiSolver->solve((Vort - FV::Div_a_Laplace_perp(1/Bsq,Pi_sum))/n_B2,
phi_boundary);
}

// Ensure that potential is set in the communication guard cells
Expand Down Expand Up @@ -464,8 +537,14 @@ void Vorticity::transform(Options &state) {
auto P = GET_NOBOUNDARY(Field3D, species["pressure"]);
auto AA = get<BoutReal>(species["AA"]);

add(species["energy_source"],
(3./2) * P * (AA / average_atomic_mass) * DivJdia);
if (boussinesq) {
add(species["energy_source"],
(3./2) * P * (AA / average_atomic_mass) * DivJdia);
} else {
auto N = GET_NOBOUNDARY(Field3D, species["density"]);
add(species["energy_source"],
(3./2) * P / N * (AA / average_atomic_mass) * DivJdia);
}
}
}

Expand Down Expand Up @@ -520,23 +599,41 @@ void Vorticity::finally(const Options &state) {
phi = get<Field3D>(state["fields"]["phi"]);

if (exb_advection) {
ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(Vort, phi, bndry_flux, poloidal_flows);

/*
if (boussinesq) {
// For now, use simplified version of ExB advection terms for Boussinesq case
ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(Vort, phi, bndry_flux, poloidal_flows);
} else {
ddt(Vort) -=
Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, bndry_flux, poloidal_flows);

// V_ExB dot Grad(Pi)
Field3D vEdotGradPi = bracket(phi, Pi, BRACKET_ARAKAWA);
vEdotGradPi.applyBoundary("free_o2");
// delp2(phi) term
Field3D DelpPhi_2B2 = 0.5 * Delp2(phi) / Bsq;
DelpPhi_2B2.applyBoundary("free_o2");

mesh->communicate(vEdotGradPi, DelpPhi_2B2);

ddt(Vort) -= FV::Div_a_Laplace_perp(0.5 / Bsq, vEdotGradPi);
*/
// V_ExB dot Grad(Pi_sum)
Field3D vEdotGradPi = bracket(phi, Pi_sum, BRACKET_ARAKAWA);
vEdotGradPi.applyBoundary("free_o2");
// delp2(phi) term
Field3D DelpPhi_2B2 = 0.5 * Delp2(phi) / Bsq;
DelpPhi_2B2.applyBoundary("free_o2");

mesh->communicate(vEdotGradPi, DelpPhi_2B2);

ddt(Vort) -= FV::Div_a_Laplace_perp(0.5 / Bsq, vEdotGradPi);

if (boussinesq) {
ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi_sum, bndry_flux,
poloidal_flows);
} else {
ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(N_sum*DelpPhi_2B2, phi, bndry_flux,
poloidal_flows);
ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, Pi_sum, bndry_flux,
poloidal_flows);

// V_ExB dot Grad(n)
Field3D vEdotGradN = bracket(phi, N_sum, BRACKET_ARAKAWA);
vEdotGradN.applyBoundary("free_o2");
mesh->communicate(vEdotGradN);

ddt(Vort) += FV::Div_a_Laplace_perp(0.5 * vEdotGradN / Bsq, phi);
}
}
}

if (state.isSection("fields") and state["fields"].isSet("DivJextra")) {
Expand Down