diff --git a/include/div_ops.hxx b/include/div_ops.hxx index 8de6aeba..096044c7 100644 --- a/include/div_ops.hxx +++ b/include/div_ops.hxx @@ -27,8 +27,8 @@ #define DIV_OPS_H #include -#include #include +#include /*! * Diffusion in index space @@ -61,9 +61,16 @@ const Field3D D4DZ4_Index(const Field3D& f); const Field2D Laplace_FV(const Field2D& k, const Field2D& f); /// Perpendicular diffusion including X and Y directions +/// Takes Div_a_Grad_perp from BOUT++ and adds flows +const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, + Field3D& flux_xlow, Field3D& flux_ylow); +/// Same but with upwinding +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f); -/// Version of function that returns flows -const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, Field3D& flux_xlow, Field3D& flux_ylow); +/// Same but with upwinding and flows +/// WARNING: Causes checkerboarding in neutral_mixed integrated test +const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, + Field3D& flux_xlow, Field3D& flux_ylow); namespace FV { @@ -105,8 +112,7 @@ struct Superbee { BoutReal sign = SIGN(gL); gL = fabs(gL); gR = fabs(gR); - BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5*gR), - BOUTMIN(gR, 0.5*gL)); + BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL)); n.L = n.c - half_slope; n.R = n.c + half_slope; } @@ -459,6 +465,235 @@ const Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } +/// Div ( a g Grad_perp(f) ) -- Perpendicular gradient-driven advection +/// +/// This version uses a slope limiter to calculate cell edge values of g in X, +/// the advects the upwind cell edge. +/// +/// 1st order upwinding is used in Y. +template +const Field3D Div_a_Grad_perp_limit(const Field3D& a, const Field3D& g, const Field3D& f) { + ASSERT2(a.getLocation() == f.getLocation()); + + Mesh* mesh = a.getMesh(); + + // Requires at least 2 communication guard cells in X, 1 in Y + ASSERT1(mesh->xstart >= 2); + ASSERT1(mesh->ystart >= 1); + + CellEdges cellboundary; + + Field3D result{zeroFrom(f)}; + + Coordinates* coord = f.getCoordinates(); + + // Flux in x + + for (int i = mesh->xstart - 1; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux from i to i+1 + + const BoutReal gradient = f(i + 1, j, k) - f(i, j, k); + + // Mid-point average boundary value of 'a' + const BoutReal aedge = 0.5 * (a(i + 1, j, k) + a(i, j, k)); + BoutReal gedge; + if (((i == mesh->xstart - 1) and mesh->firstX()) or + ((i == mesh->xend) and mesh->lastX())) { + // Mid-point average boundary value of 'g' + gedge = 0.5 * (g(i + 1, j, k) + g(i, j, k)); + } else if (gradient > 0) { + // Flux is from (i+1) to (i) + // Reconstruct `g` at left of (i+1, j, k) + + Stencil1D sg; + sg.m = g(i, j, k); + sg.c = g(i + 1, j, k); + sg.p = g(i + 2, j, k); + cellboundary(sg); // Calculate sg.R and sg.L + + gedge = sg.L; + } else { + // Flux is from (i) to (i+1) + // Reconstruct `g` at right of (i, j, k) + + Stencil1D sg; + sg.m = g(i - 1, j, k); + sg.c = g(i, j, k); + sg.p = g(i + 1, j, k); + cellboundary(sg); // Calculate sg.R and sg.L + + gedge = sg.R; + } + + // Flux across cell edge + const BoutReal fout = gradient * aedge * gedge + * (coord->J(i, j) * coord->g11(i, j) + + coord->J(i + 1, j) * coord->g11(i + 1, j)) + / (coord->dx(i, j) + coord->dx(i + 1, j)); + + result(i, j, k) += fout / (coord->dx(i, j) * coord->J(i, j)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + } + } + } + + // Y and Z fluxes require Y derivatives + + // Fields containing values along the magnetic field + Field3D fup(mesh), fdown(mesh); + Field3D aup(mesh), adown(mesh); + Field3D gup(mesh), gdown(mesh); + + // Values on this y slice (centre). + // This is needed because toFieldAligned may modify the field + Field3D ac = a; + Field3D gc = g; + Field3D fc = f; + + // Result of the Y and Z fluxes + Field3D yzresult(mesh); + yzresult.allocate(); + + if (f.hasParallelSlices() && a.hasParallelSlices() && g.hasParallelSlices()) { + // All inputs have yup and ydown + + fup = f.yup(); + fdown = f.ydown(); + + aup = a.yup(); + adown = a.ydown(); + + gup = g.yup(); + gdown = g.ydown(); + } else { + // At least one input doesn't have yup/ydown fields. + // Need to shift to/from field aligned coordinates + + fup = fdown = fc = toFieldAligned(f); + aup = adown = ac = toFieldAligned(a); + gup = gdown = gc = toFieldAligned(g); + yzresult.setDirectionY(YDirectionType::Aligned); + } + + // Y flux + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + + BoutReal coef_u = + 0.5 + * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) + + coord->g_23(i, j + 1) / SQ(coord->J(i, j + 1) * coord->Bxy(i, j + 1))); + + BoutReal coef_d = + 0.5 + * (coord->g_23(i, j) / SQ(coord->J(i, j) * coord->Bxy(i, j)) + + coord->g_23(i, j - 1) / SQ(coord->J(i, j - 1) * coord->Bxy(i, j - 1))); + + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.25 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / coord->dz(i, j); + + // Y derivative + BoutReal dfdy = 2. * (fup(i, j + 1, k) - fc(i, j, k)) + / (coord->dy(i, j + 1) + coord->dy(i, j)); + + BoutReal aedge = 0.5 * (ac(i, j, k) + aup(i, j + 1, k)); + BoutReal gedge; + if ((j == mesh->yend) and mesh->lastY(i)) { + // Midpoint boundary value + gedge = 0.5 * (gc(i, j, k) + gup(i, j + 1, k)); + } else if (dfdy > 0) { + // Flux from (j+1) to (j) + gedge = gup(i, j + 1, k); + } else { + // Flux from (j) to (j+1) + gedge = gc(i, j, k); + } + + BoutReal fout = aedge * gedge * 0.5 + * (coord->J(i, j) * coord->g23(i, j) + + coord->J(i, j + 1) * coord->g23(i, j + 1)) + * (dfdz - coef_u * dfdy); + + yzresult(i, j, k) = fout / (coord->dy(i, j) * coord->J(i, j)); + + // Calculate flux between j and j-1 + dfdz = 0.25 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) + / coord->dz(i, j); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (coord->dy(i, j) + coord->dy(i, j - 1)); + + aedge = 0.5 * (ac(i, j, k) + adown(i, j - 1, k)); + if ((j == mesh->ystart) and mesh->firstY(i)) { + gedge = 0.5 * (gc(i, j, k) + gdown(i, j - 1, k)); + } else if (dfdy > 0) { + gedge = gc(i, j, k); + } else { + gedge = gdown(i, j - 1, k); + } + + fout = aedge * gedge * 0.5 + * (coord->J(i, j) * coord->g23(i, j) + + coord->J(i, j - 1) * coord->g23(i, j - 1)) + * (dfdz - coef_d * dfdy); + + yzresult(i, j, k) -= fout / (coord->dy(i, j) * coord->J(i, j)); + } + } + } + + // Z flux + // Easier since all metrics constant in Z + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + // Coefficient in front of df/dy term + BoutReal coef = coord->g_23(i, j) + / (coord->dy(i, j + 1) + 2. * coord->dy(i, j) + coord->dy(i, j - 1)) + / SQ(coord->J(i, j) * coord->Bxy(i, j)); + + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; + + BoutReal gradient = + // df/dz + (fc(i, j, kp) - fc(i, j, k)) / coord->dz(i, j) + + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp)); + + BoutReal fout = gradient * 0.5*(ac(i,j,kp) + ac(i,j,k)) * + ((gradient > 0) ? gc(i, j, kp) : gc(i, j, k)); + + yzresult(i, j, k) += fout / coord->dz(i, j); + yzresult(i, j, kp) -= fout / coord->dz(i, j); + } + } + } + // Check if we need to transform back + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + } + + return result; +} + } // namespace FV #endif // DIV_OPS_H diff --git a/include/neutral_mixed.hxx b/include/neutral_mixed.hxx index f7a6d5a4..745feb22 100644 --- a/include/neutral_mixed.hxx +++ b/include/neutral_mixed.hxx @@ -42,26 +42,40 @@ private: BoutReal AA; ///< Atomic mass (proton = 1) Field3D Dnn; ///< Diffusion coefficient - Field3D DnnNn, DnnPn, DnnTn, DnnNVn; ///< Used for operators + Field3D DnnNn, DnnPn, DnnNVn; bool sheath_ydown, sheath_yup; BoutReal nn_floor; ///< Minimum Nn used when dividing NVn by Nn to get Vn. + BoutReal pn_floor; ///< Minimum Pn used when dividing Pn by Nn to get Tn. BoutReal flux_limit; ///< Diffusive flux limit BoutReal diffusion_limit; ///< Maximum diffusion coefficient bool neutral_viscosity; ///< include viscosity? + bool neutral_conduction; ///< Include heat conduction? bool evolve_momentum; ///< Evolve parallel momentum? + bool parallel_transport, perpendicular_transport; ///< Include transport in a particular direction? + bool precondition {true}; ///< Enable preconditioner? + bool lax_flux; ///< Use Lax flux for advection terms std::unique_ptr inv; ///< Laplacian inversion used for preconditioning Field3D density_source, pressure_source; ///< External input source Field3D Sn, Sp, Snv; ///< Particle, pressure and momentum source + Field3D sound_speed; ///< Sound speed for use with Lax flux + Field3D perp_nn_adv_src; ///< Source due to perpendicular advection operator + Field3D par_nn_adv_src; ///< Source due to parallel advection operator bool output_ddt; ///< Save time derivatives? bool diagnose; ///< Save additional diagnostics? + + // Flow diagnostics + Field3D particle_flow_xlow, particle_flow_ylow; + Field3D momentum_flow_xlow, momentum_flow_ylow; + Field3D energy_flow_xlow, energy_flow_ylow; + Field3D conduction_flow_xlow, conduction_flow_ylow; }; namespace { diff --git a/src/div_ops.cxx b/src/div_ops.cxx index 64ff8143..19e9bb33 100644 --- a/src/div_ops.cxx +++ b/src/div_ops.cxx @@ -792,7 +792,237 @@ const Field2D Laplace_FV(const Field2D &k, const Field2D &f) { return result; } +const Field3D Div_a_Grad_perp_flows(const Field3D& a, const Field3D& f, + Field3D &flow_xlow, + Field3D &flow_ylow) { + ASSERT2(a.getLocation() == f.getLocation()); + + Mesh* mesh = a.getMesh(); + + Field3D result{zeroFrom(f)}; + + Coordinates* coord = f.getCoordinates(); + + // Zero all flows + flow_xlow = 0.0; + flow_ylow = 0.0; + + // Flux in x + + int xs = mesh->xstart - 1; + int xe = mesh->xend; + + /* + if(mesh->firstX()) + xs += 1; + */ + /* + if(mesh->lastX()) + xe -= 1; + */ + + for (int i = xs; i <= xe; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux from i to i+1 + + BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + + result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); + result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_xlow(i + 1, j, k) = -1.0 * fout * coord->dy(i, j) * coord->dz(i, j); + } + } + } + + // Y and Z fluxes require Y derivatives + + // Fields containing values along the magnetic field + Field3D fup(mesh), fdown(mesh); + Field3D aup(mesh), adown(mesh); + + Field3D g23up(mesh), g23down(mesh); + Field3D g_23up(mesh), g_23down(mesh); + Field3D Jup(mesh), Jdown(mesh); + Field3D dyup(mesh), dydown(mesh); + Field3D dzup(mesh), dzdown(mesh); + Field3D Bxyup(mesh), Bxydown(mesh); + + // Values on this y slice (centre). + // This is needed because toFieldAligned may modify the field + Field3D fc = f; + Field3D ac = a; + + Field3D g23c = coord->g23; + Field3D g_23c = coord->g_23; + Field3D Jc = coord->J; + Field3D dyc = coord->dy; + Field3D dzc = coord->dz; + Field3D Bxyc = coord->Bxy; + + // Result of the Y and Z fluxes + Field3D yzresult(mesh); + yzresult.allocate(); + + if (f.hasParallelSlices() && a.hasParallelSlices()) { + // Both inputs have yup and ydown + + fup = f.yup(); + fdown = f.ydown(); + + aup = a.yup(); + adown = a.ydown(); + } else { + // At least one input doesn't have yup/ydown fields. + // Need to shift to/from field aligned coordinates + + fup = fdown = fc = toFieldAligned(f); + aup = adown = ac = toFieldAligned(a); + + yzresult.setDirectionY(YDirectionType::Aligned); + } + + if (bout::build::use_metric_3d) { + // 3D Metric, need yup/ydown fields. + // Requires previous communication of metrics + // -- should insert communication here? + if (!coord->g23.hasParallelSlices() || !coord->g_23.hasParallelSlices() + || !coord->dy.hasParallelSlices() || !coord->dz.hasParallelSlices() + || !coord->Bxy.hasParallelSlices() || !coord->J.hasParallelSlices()) { + throw BoutException("metrics have no yup/down: Maybe communicate in init?"); + } + + g23up = coord->g23.yup(); + g23down = coord->g23.ydown(); + + g_23up = coord->g_23.yup(); + g_23down = coord->g_23.ydown(); + + Jup = coord->J.yup(); + Jdown = coord->J.ydown(); + + dyup = coord->dy.yup(); + dydown = coord->dy.ydown(); + + dzup = coord->dz.yup(); + dzdown = coord->dz.ydown(); + + Bxyup = coord->Bxy.yup(); + Bxydown = coord->Bxy.ydown(); + + } else { + // No 3D metrics + // Need to shift to/from field aligned coordinates + g23up = g23down = g23c = toFieldAligned(coord->g23); + g_23up = g_23down = g_23c = toFieldAligned(coord->g_23); + Jup = Jdown = Jc = toFieldAligned(coord->J); + dyup = dydown = dyc = toFieldAligned(coord->dy); + dzup = dzdown = dzc = toFieldAligned(coord->dz); + Bxyup = Bxydown = Bxyc = toFieldAligned(coord->Bxy); + } + + // Y flux + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux between j and j+1 + int kp = (k + 1) % mesh->LocalNz; + int km = (k - 1 + mesh->LocalNz) % mesh->LocalNz; + + BoutReal coef = + 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23up(i, j + 1, k) / SQ(Jup(i, j + 1, k) * Bxyup(i, j + 1, k))); + + // Calculate Z derivative at y boundary + BoutReal dfdz = + 0.5 * (fc(i, j, kp) - fc(i, j, km) + fup(i, j + 1, kp) - fup(i, j + 1, km)) + / (dzc(i, j, k) + dzup(i, j + 1, k)); + + // Y derivative + BoutReal dfdy = + 2. * (fup(i, j + 1, k) - fc(i, j, k)) / (dyup(i, j + 1, k) + dyc(i, j, k)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + aup(i, j + 1, k)) + * (Jc(i, j, k) * g23c(i, j, k) + Jup(i, j + 1, k) * g23up(i, j + 1, k)) + * (dfdz - coef * dfdy); + + yzresult(i, j, k) = fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Calculate flux between j and j-1 + coef = + 0.5 + * (g_23c(i, j, k) / SQ(Jc(i, j, k) * Bxyc(i, j, k)) + + g_23down(i, j - 1, k) / SQ(Jdown(i, j - 1, k) * Bxydown(i, j - 1, k))); + + dfdz = 0.5 + * (fc(i, j, kp) - fc(i, j, km) + fdown(i, j - 1, kp) - fdown(i, j - 1, km)) + / (dzc(i, j, k) + dzdown(i, j - 1, k)); + + dfdy = 2. * (fc(i, j, k) - fdown(i, j - 1, k)) + / (dyc(i, j, k) + dydown(i, j - 1, k)); + + fout = 0.25 * (ac(i, j, k) + adown(i, j - 1, k)) + * (Jc(i, j, k) * g23c(i, j, k) + Jdown(i, j - 1, k) * g23down(i, j - 1, k)) + * (dfdz - coef * dfdy); + + yzresult(i, j, k) -= fout / (dyc(i, j, k) * Jc(i, j, k)); + + // Flow will be positive in the positive coordinate direction + flow_ylow(i, j, k) = -1.0 * fout * coord->dx(i, j) * coord->dz(i, j); + } + } + } + + // Z flux + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + for (int k = 0; k < mesh->LocalNz; k++) { + // Calculate flux between k and k+1 + int kp = (k + 1) % mesh->LocalNz; + + // Coefficient in front of df/dy term + BoutReal coef = g_23c(i, j, k) + / (dyup(i, j + 1, k) + 2. * dyc(i, j, k) + dydown(i, j - 1, k)) + / SQ(Jc(i, j, k) * Bxyc(i, j, k)); + + BoutReal fout = + 0.25 * (ac(i, j, k) + ac(i, j, kp)) + * (Jc(i, j, k) * coord->g33(i, j, k) + Jc(i, j, kp) * coord->g33(i, j, kp)) + * ( // df/dz + (fc(i, j, kp) - fc(i, j, k)) / dzc(i, j, k) + // - g_yz * df/dy / SQ(J*B) + - coef + * (fup(i, j + 1, k) + fup(i, j + 1, kp) - fdown(i, j - 1, k) + - fdown(i, j - 1, kp))); + + yzresult(i, j, k) += fout / (Jc(i, j, k) * dzc(i, j, k)); + yzresult(i, j, kp) -= fout / (Jc(i, j, kp) * dzc(i, j, kp)); + } + } + } + // Check if we need to transform back + if (f.hasParallelSlices() && a.hasParallelSlices()) { + result += yzresult; + } else { + result += fromFieldAligned(yzresult); + flow_ylow = fromFieldAligned(flow_ylow); + } + + return result; +} + // Div ( a Grad_perp(f) ) -- diffusion +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { ASSERT2(a.getLocation() == f.getLocation()); @@ -969,6 +1199,7 @@ const Field3D Div_a_Grad_perp_upwind(const Field3D& a, const Field3D& f) { /// ylow(i,j) /// /// +/// WARNING: Causes checkerboarding in neutral_mixed integrated test const Field3D Div_a_Grad_perp_upwind_flows(const Field3D& a, const Field3D& f, Field3D &flow_xlow, Field3D &flow_ylow) { diff --git a/src/neutral_boundary.cxx b/src/neutral_boundary.cxx index 72e625b2..548863d5 100644 --- a/src/neutral_boundary.cxx +++ b/src/neutral_boundary.cxx @@ -175,8 +175,7 @@ void NeutralBoundary::transform(Options& state) { // Outgoing neutral heat flux [W/m^2] // This is rearranged from Power for clarity - note definition of v_th. - BoutReal q = - 2 * nnsheath * tnsheath * v_th // Incident energy + BoutReal q = 2 * nnsheath * tnsheath * v_th // Incident energy - (target_energy_refl_factor * target_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy - (1 - target_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy @@ -198,11 +197,11 @@ void NeutralBoundary::transform(Options& state) { } } - // SOL edge + // SOL edge if (sol) { - if(mesh->lastX()){ // Only do this for the processor which has the edge region - for(int iy=0; iy < mesh->LocalNy ; iy++){ - for(int iz=0; iz < mesh->LocalNz; iz++){ + if (mesh->lastX()) { // Only do this for the processor which has the edge region + for (int iy = 0; iy < mesh->LocalNy; iy++) { + for (int iz = 0; iz < mesh->LocalNz; iz++) { auto i = indexAt(Nn, mesh->xend, iy, iz); // Final domain cell auto ig = indexAt(Nn, mesh->xend+1, iy, iz); // Guard cell @@ -219,11 +218,9 @@ void NeutralBoundary::transform(Options& state) { // Outgoing neutral heat flux [W/m^2] // This is rearranged from Power for clarity - note definition of v_th. - BoutReal q = - 2 * nnsheath * tnsheath * v_th // Incident energy - - (target_energy_refl_factor * target_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy - - (1 - target_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy - + BoutReal q = 2 * nnsheath * tnsheath * v_th // Incident energy + - (sol_energy_refl_factor * sol_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy + - (1 - sol_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy // Multiply by radial cell area to get power // Expanded form of the calculation for clarity @@ -245,7 +242,6 @@ void NeutralBoundary::transform(Options& state) { // Subtract from cell next to boundary energy_source[i] -= cooling_source; wall_energy_source[i] -= cooling_source; - } } } @@ -254,12 +250,12 @@ void NeutralBoundary::transform(Options& state) { // PFR edge if (pfr) { if ((mesh->firstX()) and (!mesh->periodicY(mesh->xstart))) { // do loop if inner edge and not periodic (i.e. PFR) - for(int iy=0; iy < mesh->LocalNy ; iy++){ - for(int iz=0; iz < mesh->LocalNz; iz++){ + for (int iy = 0; iy < mesh->LocalNy ; iy++) { + for (int iz = 0; iz < mesh->LocalNz; iz++) { auto i = indexAt(Nn, mesh->xstart, iy, iz); // Final domain cell auto ig = indexAt(Nn, mesh->xstart-1, iy, iz); // Guard cell - + // Calculate midpoint values at wall const BoutReal nnsheath = 0.5 * (Nn[ig] + Nn[i]); const BoutReal tnsheath = 0.5 * (Tn[ig] + Tn[i]); @@ -267,16 +263,14 @@ void NeutralBoundary::transform(Options& state) { // Thermal speed of static Maxwellian in one direction const BoutReal v_th = 0.25 * sqrt( 8*tnsheath / (PI*AA) ); // Stangeby p.69 eqns. 2.21, 2.24 - // Approach adapted from D. Power thesis 2023 + // Approach adapted from D. Power thesis 2023 BoutReal T_FC = 3 / Tnorm; // Franck-Condon temp (hardcoded for now) // Outgoing neutral heat flux [W/m^2] // This is rearranged from Power for clarity - note definition of v_th. - BoutReal q = - 2 * nnsheath * tnsheath * v_th // Incident energy - - (target_energy_refl_factor * target_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy - - (1 - target_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy - + BoutReal q = 2 * nnsheath * tnsheath * v_th // Incident energy + - (pfr_energy_refl_factor * pfr_fast_refl_fraction ) * 2 * nnsheath * tnsheath * v_th // Fast reflected energy + - (1 - pfr_fast_refl_fraction) * T_FC * nnsheath * v_th; // Thermal reflected energy // Multiply by radial cell area to get power // Expanded form of the calculation for clarity @@ -298,7 +292,6 @@ void NeutralBoundary::transform(Options& state) { // Subtract from cell next to boundary energy_source[i] -= cooling_source; wall_energy_source[i] -= cooling_source; - } } } diff --git a/src/neutral_mixed.cxx b/src/neutral_mixed.cxx index a2f848a9..226cdf0a 100644 --- a/src/neutral_mixed.cxx +++ b/src/neutral_mixed.cxx @@ -1,16 +1,18 @@ #include -#include #include #include +#include #include #include "../include/div_ops.hxx" -#include "../include/neutral_mixed.hxx" #include "../include/hermes_build_config.hxx" +#include "../include/neutral_mixed.hxx" using bout::globals::mesh; +using ParLimiter = FV::Upwind; + NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* solver) : name(name) { AUTO_TRACE(); @@ -31,16 +33,16 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* // Evolving variables e.g name is "h" or "h+" solver->add(Nn, std::string("N") + name); solver->add(Pn, std::string("P") + name); - evolve_momentum = options["evolve_momentum"] - .doc("Evolve parallel neutral momentum?") - .withDefault(true); + .doc("Evolve parallel neutral momentum?") + .withDefault(true); if (evolve_momentum) { solver->add(NVn, std::string("NV") + name); } else { - output_warn.write("WARNING: Not evolving neutral parallel momentum. NVn and Vn set to zero\n"); + output_warn.write( + "WARNING: Not evolving neutral parallel momentum. NVn and Vn set to zero\n"); NVn = 0.0; Vn = 0.0; } @@ -56,24 +58,43 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* nn_floor = options["nn_floor"] .doc("A minimum density used when dividing NVn by Nn. " "Normalised units.") - .withDefault(1e-5); + .withDefault(1e-8); + + pn_floor = nn_floor * (1./get(alloptions["units"]["eV"])); precondition = options["precondition"] .doc("Enable preconditioning in neutral model?") .withDefault(true); - flux_limit = options["flux_limit"] - .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") - .withDefault(0.2); + lax_flux = options["lax_flux"] + .doc("Enable stabilising lax flux?") + .withDefault(true); + + flux_limit = + options["flux_limit"] + .doc("Limit diffusive fluxes to fraction of thermal speed. <0 means off.") + .withDefault(0.2); diffusion_limit = options["diffusion_limit"] - .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") - .withDefault(-1.0) - / (meters * meters / seconds); // Normalise + .doc("Upper limit on diffusion coefficient [m^2/s]. <0 means off") + .withDefault(-1.0) + / (meters * meters / seconds); // Normalise neutral_viscosity = options["neutral_viscosity"] - .doc("Include neutral gas viscosity?") - .withDefault(true); + .doc("Include neutral gas viscosity?") + .withDefault(true); + + neutral_conduction = options["neutral_conduction"] + .doc("Include neutral gas heat conduction?") + .withDefault(true); + + perpendicular_transport = options["perpendicular_transport"] + .doc("Solve perpendicular transport?") + .withDefault(true); + + parallel_transport = options["parallel_transport"] + .doc("Solve parallel transport?") + .withDefault(true); if (precondition) { inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); @@ -98,10 +119,11 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* density_source = 0.0; mesh->get(density_source, std::string("N") + name + "_src"); // Allow the user to override the source - density_source = alloptions[std::string("N") + name]["source"] - .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) - .withDefault(density_source) - / (Nnorm * Omega_ci); + density_source = + alloptions[std::string("N") + name]["source"] + .doc("Source term in ddt(N" + name + std::string("). Units [m^-3/s]")) + .withDefault(density_source) + / (Nnorm * Omega_ci); // Try to read the pressure source from the mesh // Units of Pascals per second @@ -109,18 +131,22 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* mesh->get(pressure_source, std::string("P") + name + "_src"); // Allow the user to override the source pressure_source = alloptions[std::string("P") + name]["source"] - .doc(std::string("Source term in ddt(P") + name - + std::string("). Units [N/m^2/s]")) - .withDefault(pressure_source) - / (SI::qe * Nnorm * Tnorm * Omega_ci); + .doc(std::string("Source term in ddt(P") + name + + std::string("). Units [N/m^2/s]")) + .withDefault(pressure_source) + / (SI::qe * Nnorm * Tnorm * Omega_ci); // Set boundary condition defaults: Neumann for all but the diffusivity. // The dirichlet on diffusivity ensures no radial flux. // NV and V are ignored as they are hardcoded in the parallel BC code. - alloptions[std::string("Dnn") + name]["bndry_all"] = alloptions[std::string("Dnn") + name]["bndry_all"].withDefault("dirichlet"); - alloptions[std::string("T") + name]["bndry_all"] = alloptions[std::string("T") + name]["bndry_all"].withDefault("neumann"); - alloptions[std::string("P") + name]["bndry_all"] = alloptions[std::string("P") + name]["bndry_all"].withDefault("neumann"); - alloptions[std::string("N") + name]["bndry_all"] = alloptions[std::string("N") + name]["bndry_all"].withDefault("neumann"); + alloptions[std::string("Dnn") + name]["bndry_all"] = + alloptions[std::string("Dnn") + name]["bndry_all"].withDefault("dirichlet"); + alloptions[std::string("T") + name]["bndry_all"] = + alloptions[std::string("T") + name]["bndry_all"].withDefault("neumann"); + alloptions[std::string("P") + name]["bndry_all"] = + alloptions[std::string("P") + name]["bndry_all"].withDefault("neumann"); + alloptions[std::string("N") + name]["bndry_all"] = + alloptions[std::string("N") + name]["bndry_all"].withDefault("neumann"); // Pick up BCs from input file Dnn.setBoundary(std::string("Dnn") + name); @@ -134,12 +160,11 @@ NeutralMixed::NeutralMixed(const std::string& name, Options& alloptions, Solver* logPnlim.setBoundary(std::string("P") + name); Nnlim.setBoundary(std::string("N") + name); - // Product of Dnn and another parameter has same BC as Dnn - see eqns to see why this is necessary + // Product of Dnn and another parameter has same BC as Dnn - see eqns to see why this is + // necessary DnnNn.setBoundary(std::string("Dnn") + name); DnnPn.setBoundary(std::string("Dnn") + name); - DnnTn.setBoundary(std::string("Dnn") + name); DnnNVn.setBoundary(std::string("Dnn") + name); - } void NeutralMixed::transform(Options& state) { @@ -165,7 +190,7 @@ void NeutralMixed::transform(Options& state) { Vn.applyBoundary("neumann"); Vnlim.applyBoundary("neumann"); - Pnlim = floor(Nnlim * Tn, 1e-8); + Pnlim = floor(Pn, pn_floor); Pnlim.applyBoundary(); ///////////////////////////////////////////////////// @@ -264,7 +289,8 @@ void NeutralMixed::finally(const Options& state) { BoutReal neutral_lmax = 0.1 / get(state["units"]["meters"]); // Normalised length - Field3D Rnn = sqrt(Tn / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] + Field3D Rnn = + sqrt(floor(Tn, 1e-5) / AA) / neutral_lmax; // Neutral-neutral collisions [normalised frequency] if (localstate.isSet("collision_frequency")) { // Dnn = Vth^2 / sigma @@ -276,11 +302,8 @@ void NeutralMixed::finally(const Options& state) { if (flux_limit > 0.0) { // Apply flux limit to diffusion, // using the local thermal speed and pressure gradient magnitude - Field3D Dmax = flux_limit * sqrt(Tn / AA) / - (abs(Grad(logPnlim)) + 1. / neutral_lmax); - BOUT_FOR(i, Dmax.getRegion("RGN_NOBNDRY")) { - Dnn[i] = BOUTMIN(Dnn[i], Dmax[i]); - } + Field3D Dmax = flux_limit * sqrt(Tn / AA) / (abs(Grad(logPnlim)) + 1. / neutral_lmax); + BOUT_FOR(i, Dmax.getRegion("RGN_NOBNDRY")) { Dnn[i] = BOUTMIN(Dnn[i], Dmax[i]); } } if (diffusion_limit > 0.0) { @@ -295,19 +318,20 @@ void NeutralMixed::finally(const Options& state) { Dnn.applyBoundary(); // Neutral diffusion parameters have the same boundary condition as Dnn - DnnPn = Dnn * Pn; + DnnNn = Dnn * Nnlim; + DnnPn = Dnn * Pnlim; + DnnNVn = Dnn * NVn; + DnnPn.applyBoundary(); - DnnNn = Dnn * Nn; DnnNn.applyBoundary(); - Field3D DnnNVn = Dnn * NVn; DnnNVn.applyBoundary(); if (sheath_ydown) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { Dnn(r.ind, mesh->ystart - 1, jz) = -Dnn(r.ind, mesh->ystart, jz); - DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); DnnNn(r.ind, mesh->ystart - 1, jz) = -DnnNn(r.ind, mesh->ystart, jz); + DnnPn(r.ind, mesh->ystart - 1, jz) = -DnnPn(r.ind, mesh->ystart, jz); DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); } } @@ -317,22 +341,35 @@ void NeutralMixed::finally(const Options& state) { for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { Dnn(r.ind, mesh->yend + 1, jz) = -Dnn(r.ind, mesh->yend, jz); - DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); DnnNn(r.ind, mesh->yend + 1, jz) = -DnnNn(r.ind, mesh->yend, jz); + DnnPn(r.ind, mesh->yend + 1, jz) = -DnnPn(r.ind, mesh->yend, jz); DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); } } } // Sound speed appearing in Lax flux for advection terms - Field3D sound_speed = sqrt(Tn * (5. / 3) / AA); + sound_speed = 0; + if (lax_flux) { + sound_speed = sqrt(Tn * (5. / 3) / AA); + } ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = -FV::Div_par_mod(Nn, Vn, sound_speed) // Advection - + FV::Div_a_Grad_perp(DnnNn, logPnlim) // Perpendicular diffusion - ; + + perp_nn_adv_src = Div_a_Grad_perp_flows(DnnNn, logPnlim, + particle_flow_xlow, + particle_flow_ylow); // Perpendicular advection + + par_nn_adv_src = FV::Div_par_mod(Nn, Vn, sound_speed); // Parallel advection + + if (parallel_transport) { + ddt(Nn) -= par_nn_adv_src; + } + if (perpendicular_transport) { + ddt(Nn) += perp_nn_adv_src; + } Sn = density_source; // Save for possible output if (localstate.isSet("density_source")) { @@ -346,14 +383,19 @@ void NeutralMixed::finally(const Options& state) { // Neutral momentum TRACE("Neutral momentum"); - ddt(NVn) = - -AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed) // Momentum flow - - Grad_par(Pn) // Pressure gradient - + FV::Div_a_Grad_perp(DnnNVn, logPnlim) // Perpendicular diffusion - ; + if (parallel_transport) { + ddt(NVn) -= AA * FV::Div_par_fvv(Nnlim, Vn, sound_speed); // Momentum flow + ddt(NVn) -= Grad_par(Pn); // Pressure gradient + } + if (perpendicular_transport) { + ddt(NVn) += Div_a_Grad_perp_flows(DnnNVn, logPnlim, + momentum_flow_xlow, + momentum_flow_ylow); // Perpendicular advection + } + if (neutral_viscosity) { - // NOTE: The following viscosity terms are are not (yet) balanced + // NOTE: The following viscosity terms are not (yet) balanced // by a viscous heating term // Relationship between heat conduction and viscosity for neutral @@ -363,9 +405,12 @@ void NeutralMixed::finally(const Options& state) { // eta_n = (2. / 5) * kappa_n; // - ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn) // Perpendicular viscosity - + AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity - ; + if (parallel_transport) { + ddt(NVn) += AA * FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn); // Parallel viscosity + } + if (perpendicular_transport) { + ddt(NVn) += AA * FV::Div_a_Grad_perp((2. / 5) * DnnNn, Vn); // Perpendicular viscosity + } } if (localstate.isSet("momentum_source")) { @@ -382,13 +427,34 @@ void NeutralMixed::finally(const Options& state) { // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = -FV::Div_par_mod(Pn, Vn, sound_speed) // Advection - - (2. / 3) * Pn * Div_par(Vn) // Compression - + FV::Div_a_Grad_perp(DnnPn, logPnlim) // Perpendicular diffusion - + FV::Div_a_Grad_perp(DnnNn, Tn) // Conduction - + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction - ; + if (parallel_transport) { + ddt(NVn) -= FV::Div_par_mod(Pn, Vn, sound_speed); // Parallel advection + ddt(NVn) -= (2. / 3) * Pn * Div_par(Vn) ; // Parallel advection + } + if (perpendicular_transport) { + ddt(NVn) += Div_a_Grad_perp_flows(DnnPn, logPnlim, + energy_flow_xlow, energy_flow_ylow); // Perpendicular viscosity + } + + + // The factor here is 5/2 as we're advecting internal energy and pressure. + energy_flow_xlow *= 5/2; + energy_flow_ylow *= 5/2; + + if (neutral_conduction) { + if (parallel_transport) { + ddt(Pn) += FV::Div_par_K_Grad_par(DnnNn, Tn); // Parallel conduction + } + if (perpendicular_transport) { + ddt(Pn) += Div_a_Grad_perp_flows(DnnNn, Tn, + conduction_flow_xlow, conduction_flow_ylow); // Perpendicular conduction + } + } + // The factor here is likely 3/2 as this is pure energy flow, but needs checking. + conduction_flow_xlow *= 3/2; + conduction_flow_ylow *= 3/2; + Sp = pressure_source; if (localstate.isSet("energy_source")) { Sp += (2. / 3) * get(localstate["energy_source"]); @@ -396,10 +462,10 @@ void NeutralMixed::finally(const Options& state) { ddt(Pn) += Sp; BOUT_FOR(i, Pn.getRegion("RGN_ALL")) { - if ((Pn[i] < 1e-9) && (ddt(Pn)[i] < 0.0)) { + if ((Pn[i] < pn_floor * 1e-2) && (ddt(Pn)[i] < 0.0)) { ddt(Pn)[i] = 0.0; } - if ((Nn[i] < 1e-7) && (ddt(Nn)[i] < 0.0)) { + if ((Nn[i] < nn_floor * 1e-2) && (ddt(Nn)[i] < 0.0)) { ddt(Nn)[i] = 0.0; } } @@ -433,6 +499,7 @@ void NeutralMixed::outputVars(Options& state) { auto Tnorm = get(state["Tnorm"]); auto Omega_ci = get(state["Omega_ci"]); auto Cs0 = get(state["Cs0"]); + auto rho_s0 = get(state["rho_s0"]); const BoutReal Pnorm = SI::qe * Tnorm * Nnorm; state[std::string("N") + name].setAttributes({{"time_dimension", "t"}, @@ -531,7 +598,103 @@ void NeutralMixed::outputVars(Options& state) { {"long_name", name + " pressure source"}, {"species", name}, {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("S") + name + std::string("_perp_adv")], perp_nn_adv_src, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "density source due to perp advection"}, + {"long_name", name + " number density source due to perp advection"}, + {"species", name}, + {"source", "neutral_mixed"}}); + set_with_attrs(state[std::string("S") + name + std::string("_par_adv")], par_nn_adv_src, + {{"time_dimension", "t"}, + {"units", "m^-3 s^-1"}, + {"conversion", Nnorm * Omega_ci}, + {"standard_name", "density source due to par advection"}, + {"long_name", name + " number density source due to par advection"}, + {"species", name}, + {"source", "neutral_mixed"}}); + if (particle_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("ParticleFlow_") + name + std::string("_xlow")], particle_flow_xlow, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, + {"standard_name", "particle flow"}, + {"long_name", name + " particle flow in X. Note: May be incomplete."}, + {"species", name}, + {"source", "neutral_mixed"}}); + } + if (particle_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("ParticleFlow_") + name + std::string("_ylow")], particle_flow_ylow, + {{"time_dimension", "t"}, + {"units", "s^-1"}, + {"conversion", rho_s0 * SQ(rho_s0) * Nnorm * Omega_ci}, + {"standard_name", "particle flow"}, + {"long_name", name + " particle flow in Y. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_density"}}); + } + if (momentum_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("MomentumFlow_") + name + std::string("_xlow")], momentum_flow_xlow, + {{"time_dimension", "t"}, + {"units", "N"}, + {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum flow"}, + {"long_name", name + " momentum flow in X. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_momentum"}}); + } + if (momentum_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("MomentumFlow_") + name + std::string("_ylow")], momentum_flow_ylow, + {{"time_dimension", "t"}, + {"units", "N"}, + {"conversion", rho_s0 * SQ(rho_s0) * SI::Mp * Nnorm * Cs0 * Omega_ci}, + {"standard_name", "momentum flow"}, + {"long_name", name + " momentum flow in Y. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_momentum"}}); + } + if (energy_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("EnergyFlow_") + name + std::string("_xlow")], energy_flow_xlow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " power through X cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } + if (energy_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("EnergyFlow_") + name + std::string("_ylow")], energy_flow_ylow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " power through Y cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } + if (conduction_flow_xlow.isAllocated()) { + set_with_attrs(state[std::string("ConductionFlow_") + name + std::string("_xlow")],conduction_flow_xlow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " conducted power through X cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } + if (conduction_flow_ylow.isAllocated()) { + set_with_attrs(state[std::string("ConductionFlow_") + name + std::string("_ylow")], conduction_flow_ylow, + {{"time_dimension", "t"}, + {"units", "W"}, + {"conversion", rho_s0 * SQ(rho_s0) * Pnorm * Omega_ci}, + {"standard_name", "power"}, + {"long_name", name + " conducted power through Y cell face. Note: May be incomplete."}, + {"species", name}, + {"source", "evolve_pressure"}}); + } } }