diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..f8a0800 --- /dev/null +++ b/.clang-format @@ -0,0 +1,111 @@ +--- +# Note: commented out lines require clang-format >= 6 +Language: Cpp +BasedOnStyle: LLVM +AccessModifierOffset: -2 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +# AlignEscapedNewlines: Right +AlignEscapedNewlinesLeft: true +AlignOperands: true +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: false +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: false +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: true +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: false + AfterControlStatement: false + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false +# AfterExternBlock: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false +# SplitEmptyFunction: true +# SplitEmptyRecord: true +# SplitEmptyNamespace: true +BreakBeforeBinaryOperators: NonAssignment +BreakBeforeBraces: Attach +# BreakBeforeInheritanceComma: false +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +# BreakConstructorInitializers: BeforeColon +#BreakAfterJavaFieldAnnotations: false +#BreakStringLiterals: true +ColumnLimit: 90 +CommentPragmas: '^ IWYU pragma:' +# CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +# FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +# IncludeBlocks: Regroup +IncludeCategories: + - Regex: '^(<|")bout/' + Priority: 2 + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + - Regex: '^<.*>' + Priority: 4 + - Regex: '.*' + Priority: 1 +#IncludeIsMainRegex: '(Test)?$' +IndentCaseLabels: false +# IndentPPDirectives: None +IndentWidth: 2 +IndentWrappedFunctionNames: false +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +# PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 60 +PointerAlignment: Left +ReflowComments: true +SortIncludes: true +# SortUsingDeclarations: true +SpaceAfterCStyleCast: false +#SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: ControlStatements +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Cpp11 +TabWidth: 8 +UseTab: Never +... diff --git a/.gitignore b/.gitignore index 38425c0..775465a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ *.o *.nc BOUT.log.* - +*.a diff --git a/blob2d/BOUT.inp b/blob2d/BOUT.inp new file mode 100644 index 0000000..28ac6ae --- /dev/null +++ b/blob2d/BOUT.inp @@ -0,0 +1,355 @@ + +NOUT = 50 +TIMESTEP = 1e1 + +TwistShift = false # use twist-shift condition? +ballooning = false +shiftinitial = false +ShiftWithoutTwist=false + +[restart] +init_missing=true + +[mesh] +MYG=0 + +nx = 260 # Radial resolution including 4 guard cells +ny = 1 # Parallel direction +nz = 256 # number of points in toroidal direction + +Lrad = 0.3 # Radial width of domain [m] +Lpol = 0.3 # Poloidal length of domain [m] +Ltor = 100 + +Bpxy = .5 # Poloidal magnetic field [T] +Btxy = 0.0 # Toroidal magnetic field [T] +Rxy = 1.5 # Major radius [m] + + +# These indices mean whole domain in SOL (open field lines) +ixseps1 = -1 +ixseps2 = -1 + +# zShift is the toroidal angle change when following a field line +# First work out the distance moved toroidally, then convert to angle +# by dividing by major radius +#zShift = y/(2π) * Lpol * Btxy / Bpxy / Rxy + +##### + +hthe = 1.0 # This is a scale factor which can be chosen to be 1 +dy = Ltor / ny # dy * hthe is the poloidal cell size +dx = Lrad * Rxy * Bpxy / (nx - 4) # Poloidal flux +dz = Lpol/ (Rxy*nz) + +Bxy = sqrt(Bpxy^2 + Btxy^2) + +bxcvz = 1./Rxy^2 # Curvature + +paralleltransform = identity + +symmetricGlobalX = true + +shiftangle_from_zshift=false + +################################################## +# derivative methods + +[ddx] + +first = C2 +second = C2 +upwind = W3 + +[ddy] + +first = C2 +second = C2 +upwind = W3 + +[ddz] + +first = C2 +second = C2 +upwind = W3 + +################################################### +# Time-integration solver + +[solver] + +# Note: If evolving neutrals, need preconditioning +#type = cvode +#use_precon = true + +ATOL = 1.0e-10 # absolute tolerance +RTOL = 1.0e-5 # relative tolerance +mxstep = 1000000 # Maximum internal steps per output + + +[laplace] # This is used for Delp2 operator +all_terms = true +nonuniform=true + +################################################## +# Electrostatic potential solver + +[phiSolver] +inner_boundary_flags = 1 +outer_boundary_flags = 16 # INVERT_SET (2.8 * Te) + +all_terms = true +nonuniform=true # NOTE: Necessary to avoid numerical instability + +[laplacexz] +type = petsc + +rtol=1e-8 +atol=1e-12 + +ksptype = gmres # Linear iterative method + +pctype = hypre # Preconditioner. Direct "lu" or "ilu"; iterative "jacobi", "sor" + +# Set package to perform factorisation for direct solves +# "petsc" built-in solver only serial +# "superlu", "superlu_dist", "mumps", "cusparse" +factor_package = superlu_dist + +inner_boundary_flags = 16 +outer_boundary_flags = 16 + +[laplacexy] # 2D solver in X-Y +pctype = sor # Preconditioner + +atol = 1e-12 +rtol = 1e-8 + +core_bndry_dirichlet = false +pf_bndry_dirichlet = true +y_bndry_dirichlet = false + +include_y_derivs = true + +[aparSolver] +inner_boundary_flags = 0 +outer_boundary_flags = 0 + +all_terms = true +nonuniform=true + +# general settings for the model + +[Hermes] + +####################### +# Output variables + +output_ddt = false # Output time derivatives +verbose = true # Output additional fields + +####################### +# radial buffers +radial_buffers = false +radial_inner_width = 4 +radial_outer_width = 4 +radial_buffer_D = 1.0 + +resistivity_boundary = 1e-2 +resistivity_boundary_width = 0 + +####################### +# Numerical dissipation + +vepsi_dissipation = false # Parallel dissipation on Ve-Vi +vort_dissipation = false + +# Flux limiters +kappa_limit_alpha = -0.2 # SOLPS style heat flux limiter (< 0 -> off) +eta_limit_alpha = -0.5 # SOLPS style viscosity limiter (< 0 -> off) + +####################### +# Electric field and Ohm's law + +electromagnetic = false # Electromagnetic? Otherwise electrostatic +FiniteElMass = false # Finite electron mass? + +# Electrostatic potential + +newXZsolver = false # Use the XZ solver? (false -> phiSolver) +split_n0 = false # Solve n=0 separately? +split_n0_psi = false # Solve n=0 flux separately? +phi_boundary_relax = true # Relax boundaries towards zero gradient + +# NOTE: all currents switched off for fluid run +j_diamag = true # Diamagnetic current: Vort <-> Pe +j_par = false # Parallel current: Vort <-> Psi + +j_pol_terms = false # Additional terms in the vorticity equation? + +pe_par = false # Parallel pressure gradient: Pe <-> Psi +resistivity = false # Resistivity: Psi -> Pe +thermal_flux = false # Include thermal flux in Ohm's law? +thermal_force = false # Include thermal force in Ohm's law? +ion_velocity = false # Evolve parallel ion velocity? +electron_viscosity = false +ion_viscosity = false # Ion parallel viscosity +thermal_conduction = false +poloidal_flows = false # Include ExB flows in the X-Y plane? +parallel_flow = false +boussinesq = true # Use Boussinesq approximation? [always true currently] + +frecycle = 0.99 # Neutral gas recycling fraction + +carbon_fraction = 0.0 + +excitation = false # Hydrogen neutral excitation radiation + +## Settings for 2D parallel closures +sinks = false +sheath_closure = false +drift_wave = false +Lpar = 5000 # Parallel connection length [m] +sink_invlpar = 1 / Lpar # inverse parallel connection length [1/m] +alpha_dw = 0.0# 1e-3*h(0.3-x) + +# sinks = true +# sink_invlpar = 0.02 # 5m parallel connection length +# sheath_closure = true +# drift_wave = false + +####################### +# Transport coefficients +classical_diffusion = true # Collisional diffusion + +anomalous_D = -1 # Anomalous density diffusion [m^2/s] +anomalous_chi = -1 # Anomalous thermal diffusion [m^2/s] +anomalous_nu = -1 # Anomalous viscosity + +ion_neutral = 0.0 +neutral_friction = false # Friction between plasma and neutrals + +# Radial boundary fluxes +ne_bndry_flux = true +pe_bndry_flux = true +vort_bndry_flux = false + +ramp_mesh = false +ramp_timescale = 1e4 + +####################### +# Plasma sheath +sheath_model = 0 # 0 = Bohn, 1 = Loizu, 2 = Bohm + free density +sheath_yup = false +sheath_ydown = false +sheath_gamma_e = 4 # Electron sheath heat transmission +sheath_gamma_i = 2.5 # Ion sheath heat transmission +neutral_gamma = 0.0 + +startprofiles = false + +core_sources = false # Only sources in the core +adapt_source = false # Feedback on profiles (PI controller) +energy_source = true # Constant power per particle (at fixed x). False -> Constant power per volume +source_p = 1e-2 # Rate based on current error (1/timescale) +source_i = 1e-6 # Integral of error (1/timescale^2) +source_vary_g11 = false # Vary source in poloidal angle to better match radial transport + +staggered = false + +# Normalisation parameters + +Nnorm = 1e19 # Background plasma density (m^-3) +Tnorm = 50. # Electron Temperature (eV) +Bnorm = 1. # Magnetic field [T] + +[neutral] +type = none # Neutral model: none, diffusion2d, recycling, fullvelocity, mixed +viscosity = 1 # Dynamic viscosity +bulk = 0 # Bulk (volume) viscosity +conduction = 1 +neutral_gamma = 0.0 + +nn_floor = 1e-2 # Floor applied when calculating Vn = NVn / Nn +low_n_equilibriate = -1e-4 # If n < nn_floor, evolve Tn and Vn towards plasma values + +[All] +scale = 0.0 + +bndry_all = neumann_o2 +bndry_xin = neumann_o2 +bndry_xout = neumann_o2 + +[Ne] # Electron density +scale = 1 +#function = 0.1 + 0.01*sin(3*z - x) + +function = 1 + 0.27*gauss(x-0.33, 0.21/4.)*gauss((z/(2*pi))-0.5,0.21/4.) +source = 0.0 + +[Vort] +function = 0 + +[VePsi] # Ve + 0.5*beta_e*mi_me*psi +bndry_core = zerolaplace +bndry_pf = dirichlet_o2 +bndry_xout = dirichlet_o2 + +[Pe] # Electron pressure +scale = 1 +function = 1 + 1.2*(Ne:function-1) +source = Ne:source + +[Pi] +scale = 1 +function = 1 + 1.2*(Ne:function-1) +source = Ne:source + +[Ve] + +[phi] +# Radial boundaries determined by Laplacian inversion +bndry_xin = none +bndry_xout = none + +bndry_all = dirichlet_o2 + + + +[Nn] + +scale = 5e-2 +function = 1 + +[Vn] + +scale = 0 +function = 0 + +bndry_all = dirichlet_o2 + +[Vn_x] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_y] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_z] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Pn] + +scale = 1e-5 +function = 1 + +bndry_all = neumann_o2 + +[NVn] +bndry_all = dirichlet_o2 + diff --git a/diffusion2d.cxx b/diffusion2d.cxx index 9d39daa..df5b9a3 100644 --- a/diffusion2d.cxx +++ b/diffusion2d.cxx @@ -106,7 +106,7 @@ void Diffusion2D::precon(BoutReal, BoutReal gamma, BoutReal) { // Neutral gas diffusion // Solve (1 - gamma*Dnn*Delp2)^{-1} if(!inv) { - inv = Laplacian::create(); + auto inv = Laplacian::create(); // Zero value outer boundary inv->setInnerBoundaryFlags(INVERT_DC_GRAD | INVERT_AC_GRAD); diff --git a/div_ops.cxx b/div_ops.cxx index 377cc8a..74a930a 100644 --- a/div_ops.cxx +++ b/div_ops.cxx @@ -55,14 +55,14 @@ const Field3D Div_par_diffusion_index(const Field3D &f, bool bndry_flux) { continue; } BoutReal J = - 0.5 * (coord->J(i, j) + coord->J(i, j + 1)); // Jacobian at boundary + 0.5 * (coord->J(i, j, k) + coord->J.yup()(i, j + 1, k)); // Jacobian at boundary - BoutReal gradient = f(i, j + 1, k) - f(i, j, k); + BoutReal gradient = f.yup()(i, j + 1, k) - f(i, j, k); BoutReal flux = J * gradient; - result(i, j, k) += flux / coord->J(i, j); - result(i, j + 1, k) -= flux / coord->J(i, j + 1); + result(i, j, k) += flux / coord->J(i, j, k); + result.yup()(i, j + 1, k) -= flux / coord->J.yup()(i, j + 1, k); } return result; } @@ -86,10 +86,10 @@ struct Stencil1D { }; // First order upwind for testing -void Upwind(Stencil1D &n, const BoutReal h) { n.L = n.R = n.c; } +void Upwind(Stencil1D &n) { n.L = n.R = n.c; } // Fromm method -void Fromm(Stencil1D &n, const BoutReal h) { +void Fromm(Stencil1D &n) { n.L = n.c - 0.25 * (n.p - n.m); n.R = n.c + 0.25 * (n.p - n.m); } @@ -115,7 +115,7 @@ BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c)); } -void MinMod(Stencil1D &n, const BoutReal h) { +void MinMod(Stencil1D &n) { // Choose the gradient within the cell // as the minimum (smoothest) solution BoutReal slope = minmod(n.p - n.c, n.c - n.m); @@ -124,7 +124,7 @@ void MinMod(Stencil1D &n, const BoutReal h) { } // Monotonized Central limiter (Van-Leer) -void MC(Stencil1D &n, const BoutReal h) { +void MC(Stencil1D &n) { BoutReal slope = minmod(2. * (n.p - n.c), 0.5 * (n.p - n.m), 2. * (n.c - n.m)); n.L = n.c - 0.5 * slope; @@ -267,13 +267,15 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // 2) Calculate velocities on cell faces - BoutReal vU = coord->J(i, j) * (fmp - fpp) / coord->dx(i, j); // -J*df/dx - BoutReal vD = coord->J(i, j) * (fmm - fpm) / coord->dx(i, j); // -J*df/dx + BoutReal vU = 0.5 * (coord->J(i, j, k) + coord->J(i, j , kp)) * (fmp - fpp) / + coord->dx(i, j, k); // -J*df/dx + BoutReal vD = 0.5 * (coord->J(i, j, k) + coord->J(i, j , km)) * (fmm - fpm) / + coord->dx(i, j, k); // -J*df/dx - BoutReal vR = 0.5 * (coord->J(i, j) + coord->J(i + 1, j)) * (fpp - fpm) / - coord->dz; // J*df/dz - BoutReal vL = 0.5 * (coord->J(i, j) + coord->J(i - 1, j)) * (fmp - fmm) / - coord->dz; // J*df/dz + BoutReal vR = 0.5 * (coord->J(i, j, k) + coord->J(i + 1, j, k)) * (fpp - fpm) / + coord->dz(i,j,k); // J*df/dz + BoutReal vL = 0.5 * (coord->J(i, j, k) + coord->J(i - 1, j, k)) * (fmp - fmm) / + coord->dz(i,j,k); // J*df/dz // output.write("NEW: (%d,%d,%d) : (%e/%e, %e/%e)\n", i,j,k,vL,vR, // vU,vD); @@ -292,7 +294,7 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Upwind(s, mesh->dx(i,j)); // XPPM(s, mesh->dx(i,j)); // Fromm(s, coord->dx(i, j)); - MC(s, coord->dx(i, j)); + MC(s); // Right side if ((i == mesh->xend) && (mesh->lastX())) { @@ -307,18 +309,18 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Flux in from boundary flux = vR * 0.5 * (n(i + 1, j, k) + n(i, j, k)); } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } else { // Not at a boundary if (vR > 0.0) { // Flux out into next cell BoutReal flux = vR * s.R; - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); // if(i==mesh->xend) // output.write("Setting flux (%d,%d) : %e\n", @@ -342,18 +344,18 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Flux in from boundary flux = vL * 0.5 * (n(i - 1, j, k) + n(i, j, k)); } - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i - 1, j, k) += - flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } else { // Not at a boundary if (vL < 0.0) { BoutReal flux = vL * s.L; - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i - 1, j, k) += - flux / (coord->dx(i - 1, j) * coord->J(i - 1, j)); + flux / (coord->dx(i - 1, j, k) * coord->J(i - 1, j, k)); } } @@ -368,17 +370,17 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Upwind(s, coord->dz); // XPPM(s, coord->dz); // Fromm(s, coord->dz); - MC(s, coord->dz); + MC(s); if (vU > 0.0) { - BoutReal flux = vU * s.R / (coord->J(i, j) * coord->dz); - result(i, j, k) += flux; - result(i, j, kp) -= flux; + BoutReal flux = vU * s.R; + result(i, j, k) += flux / (coord->J(i, j, k) * coord->dz(i, j, k)); + result(i, j, kp) -= flux / (coord->J(i, j, kp) * coord->dz(i, j, kp)); } if (vD < 0.0) { - BoutReal flux = vD * s.L / (coord->J(i, j) * coord->dz); - result(i, j, k) -= flux; - result(i, j, km) += flux; + BoutReal flux = vD * s.L; + result(i, j, k) -= flux / (coord->J(i, j, k) * coord->dz(i, j, k)); + result(i, j, km) += flux / (coord->J(i, j, km) * coord->dz(i, j, km)); } } FV::communicateFluxes(result); @@ -426,14 +428,14 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Average dfdy to right X boundary BoutReal f_R = - 0.5 * ((coord->g11(i + 1, j) * coord->g23(i + 1, j) / - SQ(coord->Bxy(i + 1, j))) * - dfdy(i + 1, j, k) + - (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) * - dfdy(i, j, k)); + 0.5 * ((coord->g11(i + 1, j, k) * coord->g23(i + 1, j, k) / + SQ(coord->Bxy(i + 1, j, k))) * + dfdy(i + 1, j, k) + + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) * + dfdy(i, j, k)); // Advection velocity across cell face - BoutReal Vx = 0.5 * (coord->J(i + 1, j) + coord->J(i, j)) * f_R; + BoutReal Vx = 0.5 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) * f_R; // Fromm method BoutReal flux = Vx; @@ -456,9 +458,9 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, flux *= nval; } - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i + 1, j, k) -= - flux / (coord->dx(i + 1, j) * coord->J(i + 1, j)); + flux / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); } } @@ -497,13 +499,13 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, // Average dfdx to upper Y boundary BoutReal f_U = - 0.5 * ((coord->g11(i, j + 1) * coord->g23(i, j + 1) / - SQ(coord->Bxy(i, j + 1))) * - dfdx(i, j + 1, k) + - (coord->g11(i, j) * coord->g23(i, j) / SQ(coord->Bxy(i, j))) * - dfdx(i, j, k)); + 0.5 * ((coord->g11(i, j + 1, k) * coord->g23(i, j + 1, k) / + SQ(coord->Bxy(i, j + 1, k))) * + dfdx(i, j + 1, k) + + (coord->g11(i, j, k) * coord->g23(i, j, k) / SQ(coord->Bxy(i, j, k))) * + dfdx(i, j, k)); - BoutReal Vy = -0.5 * (coord->J(i, j + 1) + coord->J(i, j)) * f_U; + BoutReal Vy = -0.5 * (coord->J(i, j + 1, k) + coord->J(i, j, k)) * f_U; if (mesh->firstY(i) && !mesh->periodicY(i) && (j == mesh->ystart - 1)) { @@ -537,8 +539,8 @@ const Field3D Div_n_bxGrad_f_B_XPPM(const Field3D &n, const Field3D &f, flux *= nval; } - yresult(i, j, k) += flux / (coord->dy(i, j) * coord->J(i, j)); - yresult(i, j + 1, k) -= flux / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + yresult(i, j, k) += flux / (coord->dy(i, j, k) * coord->J(i, j, k)); + yresult(i, j + 1, k) -= flux / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); } } } @@ -588,26 +590,33 @@ const Field3D Div_Perp_Lap_FV_Index(const Field3D &as, const Field3D &fs, BoutReal gU = fs(i, j, kp) - fs(i, j, k); // Flow right - BoutReal flux = gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) * - (coord->dx(i + 1, j) + coord->dx(i, j)) * - (as(i + 1, j, k) + as(i, j, k)); + BoutReal flux = gR * 0.25 * (coord->J(i + 1, j, k) + coord->J(i, j, k)) * + (coord->dx(i + 1, j, k) + coord->dx(i, j, k)) * + (as(i + 1, j, k) + as(i, j, k)); - result(i, j, k) += flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) += flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) * - (coord->dx(i - 1, j) + coord->dx(i, j)) * - (as(i - 1, j, k) + as(i, j, k)); + flux = gL * 0.25 * (coord->J(i - 1, j, k) + coord->J(i, j, k)) * + (coord->dx(i - 1, j, k) + coord->dx(i, j, k)) * + (as(i - 1, j, k) + as(i, j, k)); - result(i, j, k) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j, k) -= flux / (coord->dx(i, j, k) * coord->J(i, j, k)); // Flow up - flux = gU * 0.5 * (as(i, j, k) + as(i, j, kp)); - result(i, j, k) += flux; + flux = gU * 0.25 * (coord->J(i, j, k) + coord->J(i, j, kp)) * + (coord->dz(i, j, k) + coord->dz(i, j, kp)) * + (as(i, j, k) + as(i, j, kp)); + + result(i, j, k) += flux / (coord->dz(i, j, k) * coord->J(i, j, k)); - flux = gD * 0.5 * (as(i, j, k) + as(i, j, km)); - result(i, j, k) -= flux; + // Flow down + flux = gD * 0.25 * (coord->J(i, j, km) + coord->J(i, j, k)) * + (coord->dz(i, j, km) + coord->dz(i, j, k)) * + (as(i, j, km) + as(i, j, k)); + + result(i, j, k) -= flux / (coord->dz(i, j, k) * coord->J(i, j, k)); } return result; @@ -628,8 +637,8 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { BoutReal d3fdx3 = (f(i + 2, j, k) - 3. * f(i + 1, j, k) + 3. * f(i, j, k) - f(i - 1, j, k)); - BoutReal flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + BoutReal flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) * + (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; if (mesh->lastX() && (i == mesh->xend)) { // Boundary @@ -644,8 +653,8 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { - (6. / 5) * f(i - 2, j, k) // f_2 ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) * + (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; } else { // No fluxes through boundary @@ -653,8 +662,8 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { } } - result(i, j, k) += flux / (coord->J(i, j) * coord->dx(i, j)); - result(i + 1, j, k) -= flux / (coord->J(i + 1, j) * coord->dx(i + 1, j)); + result(i, j, k) += flux / (coord->J(i, j, k) * coord->dx(i, j, k)); + result(i + 1, j, k) -= flux / (coord->J(i + 1, j, k) * coord->dx(i + 1, j, k)); if (j == mesh->xstart) { // Left cell boundary, no flux through boundaries @@ -670,12 +679,12 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { + (6. / 5) * f(i + 2, j, k) // f_2 ); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) * + (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; - result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); + result(i, j, k) -= flux / (coord->J(i, j, k) * coord->dx(i, j, k)); result(i - 1, j, k) += - flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + flux / (coord->J(i - 1, j, k) * coord->dx(i - 1, j, k)); } } else { @@ -683,12 +692,12 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { d3fdx3 = (f(i + 1, j, k) - 3. * f(i, j, k) + 3. * f(i - 1, j, k) - f(i - 2, j, k)); - flux = 0.25 * (coord->dx(i, j) + coord->dx(i + 1, j)) * - (coord->J(i, j) + coord->J(i + 1, j)) * d3fdx3; + flux = 0.25 * (coord->dx(i, j, k) + coord->dx(i + 1, j, k)) * + (coord->J(i, j, k) + coord->J(i + 1, j, k)) * d3fdx3; - result(i, j, k) -= flux / (coord->J(i, j) * coord->dx(i, j)); + result(i, j, k) -= flux / (coord->J(i, j, k) * coord->dx(i, j, k)); result(i - 1, j, k) += - flux / (coord->J(i - 1, j) * coord->dx(i - 1, j)); + flux / (coord->J(i - 1, j, k) * coord->dx(i - 1, j, k)); } } } @@ -697,6 +706,15 @@ const Field3D D4DX4_FV_Index(const Field3D &f, bool bndry_flux) { return result; } +const Field3D D4DZ4_Index(const Field3D &f) { + Field3D result{emptyFrom(f)}; + + BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) { + result[i] = f[i.zpp()] - 4. * f[i.zp()] + 6. * f[i] - 4. * f[i.zm()] + f[i.zmm()]; + } + return result; +} + /*! *** USED *** * X-Y diffusion * @@ -709,52 +727,56 @@ const Field2D Laplace_FV(const Field2D &k, const Field2D &f) { result.allocate(); Coordinates *coord = mesh->getCoordinates(); - + Field2D g11_2D = DC(coord->g11); + Field2D g22_2D = DC(coord->g22); + Field2D dx_2D = DC(coord->dx); + Field2D J_2D = DC(coord->J); + for (int i = mesh->xstart; i <= mesh->xend; i++) for (int j = mesh->ystart; j <= mesh->yend; j++) { // Calculate gradients on cell faces - BoutReal gR = (coord->g11(i, j) + coord->g11(i + 1, j)) * + BoutReal gR = (g11_2D(i, j) + g11_2D(i + 1, j)) * (f(i + 1, j) - f(i, j)) / - (coord->dx(i + 1, j) + coord->dx(i, j)); + (dx_2D(i + 1, j) + dx_2D(i, j)); - BoutReal gL = (coord->g11(i - 1, j) + coord->g11(i, j)) * + BoutReal gL = (g11_2D(i - 1, j) + g11_2D(i, j)) * (f(i, j) - f(i - 1, j)) / - (coord->dx(i - 1, j) + coord->dx(i, j)); + (dx_2D(i - 1, j) + dx_2D(i, j)); - BoutReal gU = (coord->g22(i, j) + coord->g22(i, j + 1)) * + BoutReal gU = (g22_2D(i, j) + g22_2D(i, j + 1)) * (f(i, j + 1) - f(i, j)) / - (coord->dy(i, j + 1) + coord->dy(i, j)); + (dx_2D(i, j + 1) + dx_2D(i, j)); - BoutReal gD = (coord->g22(i, j - 1) + coord->g22(i, j)) * + BoutReal gD = (g22_2D(i, j - 1) + g22_2D(i, j)) * (f(i, j) - f(i, j - 1)) / - (coord->dy(i, j) + coord->dy(i, j - 1)); + (dx_2D(i, j) + dx_2D(i, j - 1)); // Flow right - BoutReal flux = gR * 0.25 * (coord->J(i + 1, j) + coord->J(i, j)) * + BoutReal flux = gR * 0.25 * (J_2D(i + 1, j) + J_2D(i, j)) * (k(i + 1, j) + k(i, j)); - result(i, j) = flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j) = flux / (dx_2D(i, j) * J_2D(i, j)); // Flow left - flux = gL * 0.25 * (coord->J(i - 1, j) + coord->J(i, j)) * + flux = gL * 0.25 * (J_2D(i - 1, j) + J_2D(i, j)) * (k(i - 1, j) + k(i, j)); - result(i, j) -= flux / (coord->dx(i, j) * coord->J(i, j)); + result(i, j) -= flux / (dx_2D(i, j) * J_2D(i, j)); // Flow up - flux = gU * 0.25 * (coord->J(i, j + 1) + coord->J(i, j)) * + flux = gU * 0.25 * (J_2D(i, j + 1) + J_2D(i, j)) * (k(i, j + 1) + k(i, j)); - result(i, j) += flux / (coord->dy(i, j) * coord->J(i, j)); + result(i, j) += flux / (dx_2D(i, j) * J_2D(i, j)); // Flow down - flux = gD * 0.25 * (coord->J(i, j - 1) + coord->J(i, j)) * + flux = gD * 0.25 * (J_2D(i, j - 1) + J_2D(i, j)) * (k(i, j - 1) + k(i, j)); - result(i, j) -= flux / (coord->dy(i, j) * coord->J(i, j)); + result(i, j) -= flux / (dx_2D(i, j) * J_2D(i, j)); } return result; } diff --git a/full-velocity.cxx b/full-velocity.cxx index 3b54947..cbea1d6 100644 --- a/full-velocity.cxx +++ b/full-velocity.cxx @@ -16,8 +16,12 @@ FullVelocity::FullVelocity(Solver *solver, Mesh *mesh, Options &options) * V_x, V_y, V_z */ - coord = mesh->getCoordinates(); - + Coordinates *coord = mesh->getCoordinates(); + Field2D dx_2D = DC(coord->dx); + Field2D dy_2D = DC(coord->dy); + Field2D J_2D = DC(coord->J); + Field2D g_22_2D = DC(coord->g_22); + options.get("gamma_ratio", gamma_ratio, 5. / 3); options.get("viscosity", neutral_viscosity, 1e-2); options.get("bulk", neutral_bulk, 1e-2); @@ -75,22 +79,20 @@ FullVelocity::FullVelocity(Solver *solver, Mesh *mesh, Options &options) Txz.allocate(); Tyr.allocate(); Tyz.allocate(); - - Coordinates *coord = mesh->getCoordinates(); for (int i = 0; i < mesh->LocalNx; i++) for (int j = mesh->ystart; j <= mesh->yend; j++) { // Central differencing of coordinates BoutReal dRdtheta, dZdtheta; if (j == mesh->ystart) { - dRdtheta = (Rxy(i, j + 1) - Rxy(i, j)) / (coord->dy(i, j)); - dZdtheta = (Zxy(i, j + 1) - Zxy(i, j)) / (coord->dy(i, j)); + dRdtheta = (Rxy(i, j + 1) - Rxy(i, j)) / (dy_2D(i, j)); + dZdtheta = (Zxy(i, j + 1) - Zxy(i, j)) / (dy_2D(i, j)); } else if (j == mesh->yend) { - dRdtheta = (Rxy(i, j) - Rxy(i, j - 1)) / (coord->dy(i, j)); - dZdtheta = (Zxy(i, j) - Zxy(i, j - 1)) / (coord->dy(i, j)); + dRdtheta = (Rxy(i, j) - Rxy(i, j - 1)) / (dy_2D(i, j)); + dZdtheta = (Zxy(i, j) - Zxy(i, j - 1)) / (dy_2D(i, j)); } else { - dRdtheta = (Rxy(i, j + 1) - Rxy(i, j - 1)) / (2. * coord->dy(i, j)); - dZdtheta = (Zxy(i, j + 1) - Zxy(i, j - 1)) / (2. * coord->dy(i, j)); + dRdtheta = (Rxy(i, j + 1) - Rxy(i, j - 1)) / (2. * dy_2D(i, j)); + dZdtheta = (Zxy(i, j + 1) - Zxy(i, j - 1)) / (2. * dy_2D(i, j)); } // Match to hthe, 1/|Grad y| @@ -102,15 +104,15 @@ FullVelocity::FullVelocity(Solver *solver, Mesh *mesh, Options &options) BoutReal dRdpsi, dZdpsi; if (i == 0) { // One-sided differences - dRdpsi = (Rxy(i + 1, j) - Rxy(i, j)) / (coord->dx(i, j)); - dZdpsi = (Zxy(i + 1, j) - Zxy(i, j)) / (coord->dx(i, j)); + dRdpsi = (Rxy(i + 1, j) - Rxy(i, j)) / (dx_2D(i, j)); + dZdpsi = (Zxy(i + 1, j) - Zxy(i, j)) / (dx_2D(i, j)); } else if (i == (mesh->LocalNx - 1)) { // One-sided differences - dRdpsi = (Rxy(i, j) - Rxy(i - 1, j)) / (coord->dx(i, j)); - dZdpsi = (Zxy(i, j) - Zxy(i - 1, j)) / (coord->dx(i, j)); + dRdpsi = (Rxy(i, j) - Rxy(i - 1, j)) / (dx_2D(i, j)); + dZdpsi = (Zxy(i, j) - Zxy(i - 1, j)) / (dx_2D(i, j)); } else { - dRdpsi = (Rxy(i + 1, j) - Rxy(i - 1, j)) / (2. * coord->dx(i, j)); - dZdpsi = (Zxy(i + 1, j) - Zxy(i - 1, j)) / (2. * coord->dx(i, j)); + dRdpsi = (Rxy(i + 1, j) - Rxy(i - 1, j)) / (2. * dx_2D(i, j)); + dZdpsi = (Zxy(i + 1, j) - Zxy(i - 1, j)) / (2. * dx_2D(i, j)); } // Match to Bp, |Grad psi|. NOTE: this only works if @@ -158,6 +160,13 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, const Field3D &Ti, const Field3D &Vi) { mesh->communicate(Nn2D, Vn2D, Pn2D); + + Coordinates *coord = mesh->getCoordinates(); + Field2D dx2D = DC(coord->dx); + Field2D dy_2D = DC(coord->dy); + Field2D J_2D = DC(coord->J); + Field2D g_22_2D = DC(coord->g_22); + Field2D Bxy2D = DC(coord->Bxy); // Navier-Stokes for axisymmetric neutral gas profiles // Nn2D, Pn2D and Tn2D are unfloored @@ -181,26 +190,28 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, for (RangeIterator idwn = mesh->iterateBndryLowerY(); !idwn.isDone(); idwn.next()) { - - if (Vn2D.y(idwn.ind, mesh->ystart) < 0.0) { - // Flowing out of domain - Vn2D.y(idwn.ind, mesh->ystart - 1) = Vn2D.y(idwn.ind, mesh->ystart); - } else { - // Flowing into domain - Vn2D.y(idwn.ind, mesh->ystart - 1) = -Vn2D.y(idwn.ind, mesh->ystart); + for (int k = 0; k < mesh->LocalNz; k++) { + + if (Vn2D.y(idwn.ind, mesh->ystart, k) < 0.0) { + // Flowing out of domain + Vn2D.y(idwn.ind, mesh->ystart - 1, k) = Vn2D.y(idwn.ind, mesh->ystart, k); + } else { + // Flowing into domain + Vn2D.y(idwn.ind, mesh->ystart - 1, k) = -Vn2D.y(idwn.ind, mesh->ystart, k); + } + // Neumann boundary condition on X and Z components + Vn2D.x(idwn.ind, mesh->ystart - 1, k) = Vn2D.x(idwn.ind, mesh->ystart, k); + Vn2D.z(idwn.ind, mesh->ystart - 1, k) = Vn2D.z(idwn.ind, mesh->ystart, k); + + // Neumann conditions on density and pressure + Nn2D(idwn.ind, mesh->ystart - 1, k) = Nn2D(idwn.ind, mesh->ystart, k); + Pn2D(idwn.ind, mesh->ystart - 1, k) = Pn2D(idwn.ind, mesh->ystart, k); } - // Neumann boundary condition on X and Z components - Vn2D.x(idwn.ind, mesh->ystart - 1) = Vn2D.x(idwn.ind, mesh->ystart); - Vn2D.z(idwn.ind, mesh->ystart - 1) = Vn2D.z(idwn.ind, mesh->ystart); - - // Neumann conditions on density and pressure - Nn2D(idwn.ind, mesh->ystart - 1) = Nn2D(idwn.ind, mesh->ystart); - Pn2D(idwn.ind, mesh->ystart - 1) = Pn2D(idwn.ind, mesh->ystart); } } // Density - ddt(Nn2D) = -Div(Vn2D, Nn2D); + ddt(Nn2D) = -DC(Div(Vn2D, Nn2D)); Field2D Nn2D_floor = floor(Nn2D, 1e-2); // Velocity @@ -212,12 +223,12 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, // Convert to cylindrical coordinates for velocity // advection term. This is to avoid Christoffel symbol // terms in curvilinear geometry - Field2D vr = Txr * Vn2D.x + Tyr * Vn2D.y; // Grad R component - Field2D vz = Txz * Vn2D.x + Tyz * Vn2D.y; // Grad Z component + Field2D vr = Txr * DC(Vn2D.x) + Tyr * DC(Vn2D.y); // Grad R component + Field2D vz = Txz * DC(Vn2D.x) + Tyz * DC(Vn2D.y); // Grad Z component // Advect as scalars (no Christoffel symbols needed) - ddt(vr) = -V_dot_Grad(Vn2D, vr); - ddt(vz) = -V_dot_Grad(Vn2D, vz); + ddt(vr) = -DC(V_dot_Grad(Vn2D, vr)); + ddt(vz) = -DC(V_dot_Grad(Vn2D, vz)); // Convert back to field-aligned coordinates ddt(Vn2D).x += Urx * ddt(vr) + Uzx * ddt(vz); @@ -234,7 +245,7 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, ddt(Vn2D).x += Urx * ddt(vr) + Uzx * ddt(vz); ddt(Vn2D).y += Ury * ddt(vr) + Uzy * ddt(vz); - DivV2D = Div(Vn2D); + DivV2D = DC(Div(Vn2D)); DivV2D.applyBoundary(0.0); mesh->communicate(DivV2D); @@ -243,7 +254,7 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, ////////////////////////////////////////////////////// // Pressure - ddt(Pn2D) = -Div(Vn2D, Pn2D) - + ddt(Pn2D) = -DC(Div(Vn2D, Pn2D)) - (gamma_ratio - 1.) * Pn2D * DivV2D * floor(Nn2D, 0) / Nn2D_floor + Laplace_FV(neutral_conduction, Pn2D / Nn2D); @@ -272,14 +283,14 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, BoutReal q = neutral_gamma * Nnout * Tnout * sqrt(Tnout); // Multiply by cell area to get power BoutReal heatflux = - q * (coord->J(r.ind, mesh->ystart) + coord->J(r.ind, mesh->ystart - 1)) / - (sqrt(coord->g_22(r.ind, mesh->ystart)) + - sqrt(coord->g_22(r.ind, mesh->ystart - 1))); + q * (J_2D(r.ind, mesh->ystart) + J_2D(r.ind, mesh->ystart - 1)) / + (sqrt(g_22_2D(r.ind, mesh->ystart)) + + sqrt(g_22_2D(r.ind, mesh->ystart - 1))); // Divide by volume of cell, and multiply by 2/3 to get pressure ddt(Pn2D)(r.ind, mesh->ystart) -= (2. / 3) * heatflux / - (coord->dy(r.ind, mesh->ystart) * coord->J(r.ind, mesh->ystart)); + (dy_2D(r.ind, mesh->ystart) * J_2D(r.ind, mesh->ystart)); } for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { @@ -304,14 +315,14 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, BoutReal q = neutral_gamma * Nnout * Tnout * sqrt(Tnout); // Multiply by cell area to get power BoutReal heatflux = - q * (coord->J(r.ind, mesh->yend) + coord->J(r.ind, mesh->yend + 1)) / - (sqrt(coord->g_22(r.ind, mesh->yend)) + - sqrt(coord->g_22(r.ind, mesh->yend + 1))); + q * (J_2D(r.ind, mesh->yend) + J_2D(r.ind, mesh->yend + 1)) / + (sqrt(g_22_2D(r.ind, mesh->yend)) + + sqrt(g_22_2D(r.ind, mesh->yend + 1))); // Divide by volume of cell, and multiply by 2/3 to get pressure ddt(Pn2D)(r.ind, mesh->yend) -= (2. / 3) * heatflux / - (coord->dy(r.ind, mesh->yend) * coord->J(r.ind, mesh->yend)); + (dy_2D(r.ind, mesh->yend) * J_2D(r.ind, mesh->yend)); } // Exchange of parallel momentum. This could be done @@ -319,7 +330,7 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, // Vn2D is covariant and b = e_y / (JB) to write: // // V_{||n} = b dot V_n = Vn2D.y / (JB) - Field2D Vnpar = Vn2D.y / (coord->J * coord->Bxy); + Field2D Vnpar = DC(Vn2D.y) / (J_2D * Bxy2D); ///////////////////////////////////////////////////// // Atomic processes @@ -340,7 +351,7 @@ void FullVelocity::update(const Field3D &Ne, const Field3D &Te, ddt(Pn2D) += (2. / 3) * DC(Qi); // Momentum. Note need to turn back into covariant form - ddt(Vn2D).y += DC(F) * (coord->J * coord->Bxy) / Nn2D_floor; + ddt(Vn2D).y += DC(F) * (J_2D * Bxy2D) / Nn2D_floor; // Density evolution for (auto &i : Nn2D.getRegion("RGN_ALL")) { @@ -358,11 +369,14 @@ void FullVelocity::addPressure(int x, int y, int UNUSED(z), BoutReal dpdt) { ddt(Pn2D)(x, y) += dpdt / mesh->LocalNz; } -void FullVelocity::addMomentum(int x, int y, int UNUSED(z), BoutReal dnvdt) { +void FullVelocity::addMomentum(int x, int y, int z, BoutReal dnvdt) { // Momentum added in the direction of the magnetic field // Vn2D is covariant and b = e_y / (JB) to write: // // V_{||n} = b dot V_n = Vn2D.y / (JB) + Coordinates *coord = mesh->getCoordinates(); + BoutReal J = coord->J(x, y, z); + BoutReal Bxy = coord->Bxy(x, y, z); - ddt(Vn2D.y)(x, y) += dnvdt * (coord->J(x, y) * coord->Bxy (x, y)) / Nn2D(x, y); + ddt(Vn2D.y)(x, y, z) += dnvdt * (J * Bxy) / Nn2D(x, y); } diff --git a/hermes-2.cxx b/hermes-2.cxx index 42d2600..e7867ee 100644 --- a/hermes-2.cxx +++ b/hermes-2.cxx @@ -85,24 +85,22 @@ namespace FV { } for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes - - // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) / - (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_rp = common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) / - (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_lm = common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); - for (int k = 0; k < mesh->LocalNz; k++) { + // For right cell boundaries + BoutReal common_factor = (coord->J(i, j, k) + coord->J(i, j + 1, k)) / + (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + BoutReal flux_factor_rc = common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); + BoutReal flux_factor_rp = common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + // For left cell boundaries + common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) / + (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + BoutReal flux_factor_lc = common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); + BoutReal flux_factor_lm = common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + //////////////////////////////////////////// // Reconstruct f at the cell faces // This calculates s.R and s.L for the Right and Left @@ -228,6 +226,13 @@ const Field3D ceil(const Field3D &var, BoutReal f, REGION rgn = RGN_ALL) { // Square function for vectors Field3D SQ(const Vector3D &v) { return v * v; } +/// Modifies and returns the first argument, taking the boundary from second argument +/// This is used because unfortunately Field3D::setBoundary returns void +Field3D withBoundary(Field3D &&f, const Field3D &bndry) { + f.setBoundaryTo(bndry); + return f; +} + int Hermes::init(bool restarting) { auto& opt = Options::root(); @@ -248,11 +253,19 @@ int Hermes::init(bool restarting) { j_diamag = optsc["j_diamag"] .doc("Diamagnetic current: Vort <-> Pe") .withDefault(true); + + if (optsc["j_diamag_scale"].isSet()) { + j_diamag_scale_generator = FieldFactory::get()->parse("hermes:j_diamag_scale", Options::getRoot()); + SAVE_REPEAT(j_diamag_scale); + } else { + j_diamag_scale_generator = FieldFactory::get()->parse("1"); + } j_par = optsc["j_par"] .doc("Parallel current: Vort <-> Psi") .withDefault(true); - + + OPTION(optsc, j_pol_pi, false); OPTION(optsc, parallel_flow, true); OPTION(optsc, parallel_flow_p_term, parallel_flow); OPTION(optsc, pe_par, true); @@ -265,7 +278,8 @@ int Hermes::init(bool restarting) { .withDefault(true); OPTION(optsc, electron_viscosity, true); - OPTION(optsc, ion_viscosity, true); + ion_viscosity = optsc["ion_viscosity"].doc("Include ion viscosity?").withDefault(true); + ion_viscosity_par = optsc["ion_viscosity_par"].doc("Include parallel diffusion of ion momentum?").withDefault(ion_viscosity); electron_neutral = optsc["electron_neutral"] .doc("Include electron-neutral collisions in resistivity?") @@ -293,7 +307,9 @@ int Hermes::init(bool restarting) { OPTION(optsc, pe_bndry_flux, true); OPTION(optsc, vort_bndry_flux, false); - OPTION(optsc, ramp_mesh, true); + ramp_mesh = optsc["ramp_mesh"] + .doc("Add profiles from mesh file over a period of time") + .withDefault(false); OPTION(optsc, ramp_timescale, 1e4); OPTION(optsc, energy_source, false); @@ -302,8 +318,6 @@ int Hermes::init(bool restarting) { .doc("A fixed ion-neutral collision rate, normalised to ion cyclotron frequency.") .withDefault(0.0); - OPTION(optsc, staggered, false); - OPTION(optsc, boussinesq, false); OPTION(optsc, sinks, false); @@ -337,6 +351,9 @@ int Hermes::init(bool restarting) { OPTION(optsc, floor_num_cs, -1.0); OPTION(optsc, vepsi_dissipation, false); OPTION(optsc, vort_dissipation, false); + phi_dissipation = optsc["phi_dissipation"] + .doc("Add a dissipation term to vorticity, depending on reconstruction of potential?") + .withDefault(false); OPTION(optsc, ne_hyper_z, -1.0); OPTION(optsc, pe_hyper_z, -1.0); @@ -355,6 +372,8 @@ int Hermes::init(bool restarting) { OPTION(optsc, sheath_yup, true); // Apply sheath at yup? OPTION(optsc, sheath_ydown, true); // Apply sheath at ydown? OPTION(optsc, test_boundaries, false); // Test boundary conditions + + nesheath_floor = optsc["nesheath_floor"].doc("Ne sheath lower limit").withDefault(1e-5); // Fix profiles in SOL OPTION(optsc, sol_fix_profiles, false); @@ -363,10 +382,31 @@ int Hermes::init(bool restarting) { sol_te = FieldFactory::get()->parse("sol_te", &optsc); } - OPTION(optsc, radial_buffers, false); + radial_buffers = optsc["radial_buffers"] + .doc("Turn on radial buffer regions?").withDefault(false); OPTION(optsc, radial_inner_width, 4); OPTION(optsc, radial_outer_width, 4); OPTION(optsc, radial_buffer_D, 1.0); + + // Only average in Y if in a closed field line region + radial_inner_averagey = mesh->periodicY(1) + & optsc["radial_core_averagey"] + .doc("Average Ne, Pe and Pi in Y in core radial buffer?") + .withDefault(true); + // Note: It is probably a bad idea to do this, but you never know... + radial_inner_averagey_vort = mesh->periodicY(1) + & optsc["radial_core_averagey_vort"] + .doc("Average Vort in Y in core radial buffer?") + .withDefault(false); + // These treatments of NVi might be reasonable choices, but are off by default + radial_inner_averagey_nvi = mesh->periodicY(1) + & optsc["radial_core_averagey_nvi"] + .doc("Average NVi in Y in core radial buffer?") + .withDefault(false); + radial_inner_zero_nvi = mesh->periodicY(1) + & optsc["radial_core_zero_nvi"] + .doc("Damp NVi toward zero in core radial buffer?") + .withDefault(false); resistivity_boundary = optsc["resistivity_boundary"] .doc("Normalised resistivity in radial boundary region") @@ -479,6 +519,7 @@ int Hermes::init(bool restarting) { for (int y = mesh->ystart; y <= mesh->yend; y++) { Sn(x, y) = 0.0; Spe(x, y) = 0.0; + Spi(x, y) = 0.0; } } } @@ -545,23 +586,39 @@ int Hermes::init(bool restarting) { } } - OPTION(optsc, adapt_source, false); - if (adapt_source) { - // Adaptive sources to match profiles + bool adapt_source = optsc["adapt_source"] + .doc("Vary sources in time to match core profiles. This sets the default for adapt_source_p and adapt_source_n").withDefault(false); + adapt_source_p = optsc["adapt_source_p"] + .doc("Vary pressure source in time to match core profiles").withDefault(adapt_source); + adapt_source_n = optsc["adapt_source_n"] + .doc("Vary density source in time to match core profiles").withDefault(adapt_source); + sources_positive = optsc["sources_positive"].doc("Ensure that sources are always positive").withDefault(true); + + if (adapt_source_p) { + // Adaptive pressure sources to match profiles // PI controller, including an integrated difference term OPTION(optsc, source_p, 1e-2); OPTION(optsc, source_i, 1e-6); - Field2D Snsave = copy(Sn); Field2D Spesave = copy(Spe); Field2D Spisave = copy(Spi); - SOLVE_FOR(Sn, Spe, Spi); - Sn = Snsave; + SOLVE_FOR(Spe, Spi); Spe = Spesave; Spi = Spisave; } else { - SAVE_ONCE(Sn, Spe, Spi); + SAVE_ONCE(Spe, Spi); + } + + if (adapt_source_n) { + // Adaptive density sources to match profile + OPTION(optsc, source_p, 1e-2); + OPTION(optsc, source_i, 1e-6); + Field2D Snsave = copy(Sn); + SOLVE_FOR(Sn); + Sn = Snsave; + } else { + SAVE_ONCE(Sn); } ///////////////////////////////////////////////////////// @@ -766,10 +823,24 @@ int Hermes::init(bool restarting) { } } - if (Options::root()["mesh"]["paralleltransform"].as() == "shifted") { - Field2D I; - mesh->get(I, "sinty"); - Curlb_B.z += I * Curlb_B.x; + if (Options::root()["mesh"]["paralleltransform:type"].as() == "shifted") { + // Check if the gridfile was created for "shiftedmetric" or for "identity" parallel + // transform + std::string gridfile_parallel_transform; + if (mesh->get(gridfile_parallel_transform, "parallel_transform")) { + // Did not find in gridfile, indicates older gridfile, which generated output for + // field-aligned coordinates, i.e. "identity" parallel transform + gridfile_parallel_transform = "identity"; + } + if (gridfile_parallel_transform == "identity") { + Field2D I; + mesh->get(I, "sinty"); + Curlb_B.z += I * Curlb_B.x; + } else if ((gridfile_parallel_transform != "shifted") and + (gridfile_parallel_transform != "shiftedmetric")){ + throw BoutException("Gridfile generated for unsupported parallel transform %s", + gridfile_parallel_transform.c_str()); + } } Curlb_B.x /= Bnorm; @@ -778,60 +849,100 @@ int Hermes::init(bool restarting) { Curlb_B *= 2. / coord->Bxy; - SAVE_REPEAT(phi); + ////////////////////////////////////////////////////////////// + // Electromagnetic fields + + if (j_par | j_diamag) { + // Only needed if there are any currents + SAVE_REPEAT(phi); - if (j_par) { - SAVE_REPEAT(Ve); + if (j_par) { + SAVE_REPEAT(Ve); - if (electromagnetic) - SAVE_REPEAT(psi); - } + if (electromagnetic) + SAVE_REPEAT(psi); + } - OPTION(optsc, split_n0, false); // Split into n=0 and n~=0 - OPTION(optsc, split_n0_psi, split_n0); - // Phi solver - if (phi3d) { + OPTION(optsc, split_n0, false); // Split into n=0 and n~=0 + OPTION(optsc, split_n0_psi, split_n0); + // Phi solver + if (phi3d) { #ifdef PHISOLVER - phiSolver3D = Laplace3D::create(); + phiSolver3D = Laplace3D::create(); #endif - } else { - if (split_n0) { - // Create an XY solver for n=0 component - laplacexy = new LaplaceXY(mesh); - // Set coefficients for Boussinesq solve - laplacexy->setCoefs(1. / SQ(coord->Bxy), 0.0); - phi2D = 0.0; // Starting guess - } - - // Create an XZ solver - OPTION(optsc, newXZsolver, false); - if (newXZsolver) { - // Test new LaplaceXZ solver - newSolver = LaplaceXZ::create(mesh); - // Set coefficients for Boussinesq solve - newSolver->setCoefs(1. / SQ(coord->Bxy), 0.0); } else { - // Use older Laplacian solver - phiSolver = Laplacian::create(&opt["phiSolver"]); - // Set coefficients for Boussinesq solve - phiSolver->setCoefC(1. / SQ(coord->Bxy)); + if (split_n0) { + // Create an XY solver for n=0 component + laplacexy = new LaplaceXY(mesh); + // Set coefficients for Boussinesq solve + laplacexy->setCoefs(1. / SQ(DC(coord->Bxy)), 0.0); + phi2D = 0.0; // Starting guess + } + + // Create an XZ solver + OPTION(optsc, newXZsolver, false); + if (newXZsolver) { + // Test new LaplaceXZ solver + newSolver = LaplaceXZ::create(mesh); + // Set coefficients for Boussinesq solve + newSolver->setCoefs(1. / SQ(coord->Bxy), 0.0); + } else { + // Use older Laplacian solver + phiSolver = Laplacian::create(&opt["phiSolver"]); + // Set coefficients for Boussinesq solve + phiSolver->setCoefC(1. / SQ(coord->Bxy)); + } } phi = 0.0; - } + phi.setBoundary("phi"); // For y boundaries + + phi_boundary_relax = optsc["phi_boundary_relax"] + .doc("Relax x boundaries of phi towards Neumann?") + .withDefault(false); - // Apar (Psi) solver - aparSolver = LaplaceXZ::create(mesh, &opt["aparSolver"], CELL_CENTRE); - if (split_n0_psi) { - // Use another XY solver for n=0 psi component - aparXY = new LaplaceXY(mesh); - psi2D = 0.0; - } + // Add phi to restart files so that the value in the boundaries + // is restored on restart. This is done even when phi is not evolving, + // so that phi can be saved and re-loaded + restart.addOnce(phi, "phi"); + + if (phi_boundary_relax) { + + if (!restarting) { + // Start by setting to the sheath current = 0 boundary value + + Field3D Nelim = floor(Ne, 1e-5); + Telim = floor(Pe / Nelim, 0.1 / Tnorm); + Tilim = floor(Pi / Nelim, 0.1 / Tnorm); + + phi.setBoundaryTo(DC( + (log(0.5 * sqrt(mi_me / PI)) + log(sqrt(Telim / (Telim + Tilim)))) * Telim)); + } + + // Set the last update time to -1, so it will reset + // the first time RHS function is called + phi_boundary_last_update = -1.; + + phi_boundary_timescale = optsc["phi_boundary_timescale"] + .doc("Timescale for phi boundary relaxation [seconds]") + .withDefault(1e-4) + * Omega_ci; // Normalise to internal time units + } + + // Apar (Psi) solver + aparSolver = LaplaceXZ::create(mesh, &opt["aparSolver"], CELL_CENTRE); + if (split_n0_psi) { + // Use another XY solver for n=0 psi component + aparXY = new LaplaceXY(mesh); + psi2D = 0.0; + } - Ve.setBoundary("Ve"); - nu.setBoundary("nu"); - phi.setBoundary("phi"); // For y boundaries - Jpar.setBoundary("Jpar"); + Ve.setBoundary("Ve"); + nu.setBoundary("nu"); + Jpar.setBoundary("Jpar"); + psi = 0.0; + } + nu = 0.0; kappa_epar = 0.0; Dn = 0.0; @@ -844,8 +955,10 @@ int Hermes::init(bool restarting) { SAVE_REPEAT(tau_e, tau_i); - SAVE_REPEAT(kappa_epar); // Parallel electron heat conductivity - SAVE_REPEAT(kappa_ipar); // Parallel ion heat conductivity + if (thermal_conduction || sinks) { + SAVE_REPEAT(kappa_epar); // Parallel electron heat conductivity + SAVE_REPEAT(kappa_ipar); // Parallel ion heat conductivity + } if (resistivity) { SAVE_REPEAT(nu); // Parallel resistivity @@ -861,9 +974,7 @@ int Hermes::init(bool restarting) { // Sources added to Ne, Pe and Pi equations SAVE_REPEAT(NeSource, PeSource, PiSource); } - - psi = phi = 0.0; - + // Preconditioner setPrecon((preconfunc)&Hermes::precon); @@ -883,6 +994,9 @@ int Hermes::rhs(BoutReal t) { sheath_model = 0; } + // Factor scaling the diamagnetic current + j_diamag_scale = j_diamag_scale_generator->generate(0, 0, 0, t); + // Communicate evolving variables // Note: Parallel slices are not calculated because parallel derivatives // are calculated using field aligned quantities @@ -977,7 +1091,24 @@ int Hermes::rhs(BoutReal t) { sound_speed = floor(sound_speed, floor_num_cs); } sound_speed.applyBoundary("neumann"); - + + // Maximum wave speed either electron sound speed or Alfven speed + Field3D max_speed = Bnorm * coord->Bxy / + sqrt(SI::mu0 * AA * SI::Mp * Nnorm * Nelim) / + Cs0; // Alfven speed (normalised by Cs0) + Field3D elec_sound = sqrt(mi_me) * sound_speed; // Electron sound speed + for (auto& i : max_speed.getRegion(RGN_ALL)) { + if (elec_sound[i] > max_speed[i]) { + max_speed[i] = elec_sound[i]; + } + + // Limit to 100x reference sound speed or light speed + BoutReal lim = BOUTMIN(100., 3e8/Cs0); + if (max_speed[i] > lim) { + max_speed[i] = lim; + } + } + ////////////////////////////////////////////////////////////// // Calculate electrostatic potential phi // @@ -990,6 +1121,96 @@ int Hermes::rhs(BoutReal t) { // phi = 0.0; // Already set in initialisation } else { // Solve phi from Vorticity + + // 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. + // Note: For now the boundary values are all at the midpoint, + // and only phi is considered, not phi + Pi which is handled in Boussinesq solves + Field2D phi_boundary2d; + Field3D phi_boundary3d; + + if (phi_boundary_relax) { + // Update the boundary regions by relaxing towards zero gradient + // on a given timescale. + + if (phi_boundary_last_update < 0.0) { + // First time this has been called. + phi_boundary_last_update = t; + + } else if (t > phi_boundary_last_update) { + // Only update if time has advanced + // Uses an exponential decay of the weighting of the value in the boundary + // so that the solution is well behaved for arbitrary steps + BoutReal weight = exp(-(t - phi_boundary_last_update) / phi_boundary_timescale); + phi_boundary_last_update = t; + + if (mesh->firstX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal phivalue = 0.0; + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xstart, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + + // Old value of phi at boundary + BoutReal oldvalue = + 0.5 * (phi(mesh->xstart - 1, j, 0) + phi(mesh->xstart, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = + weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xstart - 1, j, k) = 2.*newvalue - phi(mesh->xstart, j, k); + + // Note: This seems to make a difference, but don't know why. + // Without this, get convergence failures with no apparent instability + // (all fields apparently smooth, well behaved) + phi(mesh->xstart - 2, j, k) = phi(mesh->xstart - 1, j, k); + } + } + } + + if (mesh->lastX()) { + for (int j = mesh->ystart; j <= mesh->yend; j++) { + BoutReal phivalue = 0.0; + for (int k = 0; k < mesh->LocalNz; k++) { + phivalue += phi(mesh->xend, j, k); + } + phivalue /= mesh->LocalNz; // Average in Z of point next to boundary + + // Old value of phi at boundary + BoutReal oldvalue = 0.5 * (phi(mesh->xend + 1, j, 0) + phi(mesh->xend, j, 0)); + + // New value of phi at boundary, relaxing towards phivalue + BoutReal newvalue = weight * oldvalue + (1. - weight) * phivalue; + + // Set phi at the boundary to this value + for (int k = 0; k < mesh->LocalNz; k++) { + phi(mesh->xend + 1, j, k) = 2.*newvalue - phi(mesh->xend, j, k); + + // Note: This seems to make a difference, but don't know why. + // Without this, get convergence failures with no apparent instability + // (all fields apparently smooth, well behaved) + phi(mesh->xend + 2, j, k) = phi(mesh->xend + 1, j, k); + } + } + } + } + phi_boundary3d = phi; + } else { + // phi_boundary_relax = false + // + // Set boundary from temperature, to be consistent with j=0 at sheath + + // Sheath multiplier Te -> phi (2.84522 for Deuterium if Ti = 0) + phi_boundary2d = + DC((log(0.5 * sqrt(mi_me / PI)) + log(sqrt(Telim / (Telim + Tilim)))) * Telim); + + phi_boundary3d = phi_boundary2d; + } + if (phi3d) { #ifdef PHISOLVER phiSolver3D->setCoefC(Ne / SQ(coord->Bxy)); @@ -1007,33 +1228,71 @@ int Hermes::rhs(BoutReal t) { // Phi flags should be set in BOUT.inp // phiSolver->setInnerBoundaryFlags(INVERT_DC_GRAD); // phiSolver->setOuterBoundaryFlags(INVERT_SET); + + if (boussinesq) { - // Sheath multiplier Te -> phi (2.84522 for Deuterium) - BoutReal sheathmult = log(0.5 * sqrt(mi_me / PI)); + // 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) { + 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_boundary3d(mesh->xstart - 1, j, k) = + 0.5 + * (phi_boundary3d(mesh->xstart - 1, j, k) + + phi_boundary3d(mesh->xstart, j, k) + + Pi(mesh->xstart - 1, j, k) + + 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_boundary3d(mesh->xend + 1, j, k) = + 0.5 + * (phi_boundary3d(mesh->xend + 1, j, k) + + phi_boundary3d(mesh->xend, j, k) + + Pi(mesh->xend + 1, j, k) + + Pi(mesh->xend, j, k)); + } + } + } if (split_n0) { //////////////////////////////////////////// // Boussinesq, split // Split into axisymmetric and non-axisymmetric components Field2D Vort2D = DC(Vort); // n=0 component - // Set the boundary to 2.8*Te - // phi2D.setBoundaryTo(sheathmult * floor(Te.DC(), 0.0)); - phi2D.setBoundaryTo(sheathmult * DC(Telim)); + if (!phi_boundary2d.isAllocated()) { + // Make sure that the 2D boundary field is set + phi_boundary2d = DC(phi_boundary3d); + } + + // Set the boundary + phi2D.setBoundaryTo(phi_boundary2d); phi2D = laplacexy->solve(Vort2D, phi2D); // Solve non-axisymmetric part using X-Z solver if (newXZsolver) { newSolver->setCoefs(1. / SQ(coord->Bxy), 0.0); - phi = newSolver->solve(Vort - Vort2D, phi); + phi = newSolver->solve(Vort - Vort2D, + // Second argument is initial guess. Use current phi, and update boundary + withBoundary(phi + Pi - phi2D, // Value in domain + phi_boundary3d - phi_boundary2d)); // boundary } else { phiSolver->setCoefC(1. / SQ(coord->Bxy)); - // phi = phiSolver->solve((Vort-Vort2D)*SQ(coord->Bxy), phi); phi = phiSolver->solve((Vort - Vort2D) * SQ(coord->Bxy), - sheathmult * (Telim - DC(Telim))); + phi_boundary3d - phi_boundary2d); // Note: non-zero due to Pi variation } phi += phi2D; // Add axisymmetric part } else { @@ -1042,14 +1301,15 @@ int Hermes::rhs(BoutReal t) { // Solve all components using X-Z solver if (newXZsolver) { - // Use the new LaplaceXY solver - // newSolver->setCoefs(1./SQ(coord->Bxy), 0.0); // Set when - // initialised - phi = newSolver->solve(Vort, phi); + // Use the new LaplaceXZ solver + // newSolver->setCoefs(1./SQ(coord->Bxy), 0.0); // Set when initialised + phi = newSolver->solve(Vort, + withBoundary(phi + Pi, // Value in domain + phi_boundary3d)); // boundary } else { // Use older Laplacian solver // phiSolver->setCoefC(1./SQ(coord->Bxy)); // Set when initialised - phi = phiSolver->solve(Vort * SQ(coord->Bxy), (sheathmult + log(sqrt(Telim / (Telim + Tilim)))) * Telim ); + phi = phiSolver->solve(Vort * SQ(coord->Bxy), phi_boundary3d); } } @@ -1061,9 +1321,6 @@ int Hermes::rhs(BoutReal t) { // Non-Boussinesq // throw BoutException("Non-Boussinesq not implemented yet"); - - phiSolver->setCoefC(Nelim / SQ(coord->Bxy)); - phi = phiSolver->solve(Vort * SQ(coord->Bxy) / Nelim, sheathmult * Telim); } } phi.applyBoundary(t); @@ -1210,7 +1467,7 @@ int Hermes::rhs(BoutReal t) { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Zero-gradient density - BoutReal nesheath = floor(Ne(r.ind, mesh->ystart, jz), 0.0); + BoutReal nesheath = floor(Ne(r.ind, mesh->ystart, jz), nesheath_floor); // Temperature at the sheath entrance BoutReal tesheath = floor(Te(r.ind, mesh->ystart, jz), 0.0); @@ -1405,8 +1662,8 @@ int Hermes::rhs(BoutReal t) { // Zero-gradient density BoutReal nesheath = 0.5 * (3. * Ne(r.ind, mesh->ystart, jz) - Ne(r.ind, mesh->ystart + 1, jz)); - if (nesheath < 0.0) - nesheath = 0.0; + if (nesheath < nesheath_floor) + nesheath = nesheath_floor; // Temperature at the sheath entrance BoutReal tesheath = floor(Te(r.ind, mesh->ystart, jz), 0.0); @@ -1431,7 +1688,7 @@ int Hermes::rhs(BoutReal t) { -sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); // J = n*(Vi - Ve) BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-10) { + if (nesheath <= nesheath_floor) { vesheath = visheath; jsheath = 0.0; } @@ -1517,6 +1774,75 @@ int Hermes::rhs(BoutReal t) { } break; } + case 4: { + // Bohm sheath, with free boundary on the density and pressure + // The extrapolation is done using ratios i.e. linear in log(Ne) and log(P) + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + + // Free boundary, linear extrapolation of logarithms + for (int jy = mesh->ystart - 1; jy >= 0; jy--) { + Ne(r.ind, jy, jz) = SQ(Ne(r.ind, mesh->ystart, jz)) / Ne(r.ind, mesh->ystart + 1, jz); + Pe(r.ind, jy, jz) = SQ(Pe(r.ind, mesh->ystart, jz)) / Pe(r.ind, mesh->ystart + 1, jz); + Pi(r.ind, jy, jz) = SQ(Pi(r.ind, mesh->ystart, jz)) / Pi(r.ind, mesh->ystart + 1, jz); + } + + // Value at sheath from linear interpolation + // consistent with finite different methods + BoutReal nesheath = 0.5 * (Ne(r.ind, mesh->ystart, jz) + + Ne(r.ind, mesh->ystart - 1, jz)); + BoutReal tesheath = 0.5 * (Te(r.ind, mesh->ystart, jz) + + Pe(r.ind, mesh->ystart-1, jz) / Ne(r.ind, mesh->ystart-1, jz)); + BoutReal tisheath = 0.5 * (Ti(r.ind, mesh->ystart, jz) + + Pi(r.ind, mesh->ystart-1, jz) / Ne(r.ind, mesh->ystart-1, jz)); + + // Zero-gradient potential + BoutReal phisheath = phi(r.ind, mesh->ystart, jz); + + // Ion velocity goes to the sound speed + BoutReal visheath = + -sqrt(tesheath + tisheath); // Sound speed outwards + + if (Vi(r.ind, mesh->ystart, jz) < visheath) { + // If plasma is faster, go to plasma velocity + visheath = Vi(r.ind, mesh->ystart, jz); + } + + // Sheath current + BoutReal phi_te = floor(phisheath / tesheath, 0.0); + BoutReal vesheath = + -sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); + // J = n*(Vi - Ve) + BoutReal jsheath = nesheath * (visheath - vesheath); + if (nesheath < 1e-10) { + vesheath = visheath; + jsheath = 0.0; + } + + // Apply boundary condition half-way between cells + for (int jy = mesh->ystart - 1; jy >= 0; jy--) { + // Neumann conditions + phi(r.ind, jy, jz) = phisheath; + Vort(r.ind, jy, jz) = Vort(r.ind, mesh->ystart, jz); + + // Here zero-gradient Te, heat flux applied later + Te(r.ind, jy, jz) = Te(r.ind, mesh->ystart, jz); + Ti(r.ind, jy, jz) = Ti(r.ind, mesh->ystart, jz); + + // Dirichlet conditions to set flows + Vi(r.ind, jy, jz) = 2. * visheath - Vi(r.ind, mesh->ystart, jz); + Ve(r.ind, jy, jz) = 2. * vesheath - Ve(r.ind, mesh->ystart, jz); + Jpar(r.ind, jy, jz) = 2. * jsheath - Jpar(r.ind, mesh->ystart, jz); + NVi(r.ind, jy, jz) = + 2. * nesheath * visheath - NVi(r.ind, mesh->ystart, jz); + } + } + } + break; + } + default: { + throw BoutException("Sheath model %d not implemented", sheath_model); + } } } else { // sheath_ydown == false @@ -1536,6 +1862,8 @@ int Hermes::rhs(BoutReal t) { Pi(r.ind, jy, jz) = Ne(r.ind, jy, jz) * Ti(r.ind, jy, jz); Pilim(r.ind, jy, jz) = Nelim(r.ind, jy, jz) * Tilim(r.ind, jy, jz); + phi(r.ind, jy, jz) = phi(r.ind, mesh->ystart, jz); + // No flows Vi(r.ind, jy, jz) = -Vi(r.ind, mesh->ystart, jz); NVi(r.ind, jy, jz) = -NVi(r.ind, mesh->ystart, jz); @@ -1553,7 +1881,7 @@ int Hermes::rhs(BoutReal t) { for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Zero-gradient density - BoutReal nesheath = floor(Ne(r.ind, mesh->yend, jz), 0.0); + BoutReal nesheath = floor(Ne(r.ind, mesh->yend, jz), nesheath_floor); // Temperature at the sheath entrance BoutReal tesheath = floor(Te(r.ind, mesh->yend, jz), 0.0); @@ -1742,8 +2070,8 @@ int Hermes::rhs(BoutReal t) { // Zero-gradient density BoutReal nesheath = 0.5 * (3. * Ne(r.ind, mesh->yend, jz) - Ne(r.ind, mesh->yend - 1, jz)); - if (nesheath < 0.0) - nesheath = 0.0; + if (nesheath < nesheath_floor) + nesheath = nesheath_floor; // Temperature at the sheath entrance BoutReal tesheath = floor(Te(r.ind, mesh->yend, jz), 0.0); @@ -1767,7 +2095,7 @@ int Hermes::rhs(BoutReal t) { sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); // J = n*(Vi - Ve) BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-10) { + if (nesheath <= nesheath_floor) { vesheath = visheath; jsheath = 0.0; } @@ -1779,7 +2107,7 @@ int Hermes::rhs(BoutReal t) { Vort(r.ind, jy, jz) = Vort(r.ind, mesh->yend, jz); // Here zero-gradient Te, heat flux applied later - Te(r.ind, jy, jz) = tesheath; + Te(r.ind, jy, jz) = Te(r.ind, mesh->yend, jz); Ti(r.ind, jy, jz) = Ti(r.ind, mesh->yend, jz); // Dirichlet conditions @@ -1851,6 +2179,73 @@ int Hermes::rhs(BoutReal t) { } break; } + case 4: { + // Bohm sheath, with free boundary on the density and pressure + // The extrapolation is done using ratios i.e. linear in log(Ne) and log(P) + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + // Free boundary, linear extrapolation of logarithms + for (int jy = mesh->yend + 1; jy < mesh->LocalNy; jy++) { + Ne(r.ind, jy, jz) = SQ(Ne(r.ind, mesh->yend, jz)) / Ne(r.ind, mesh->yend - 1, jz); + Pe(r.ind, jy, jz) = SQ(Pe(r.ind, mesh->yend, jz)) / Pe(r.ind, mesh->yend - 1, jz); + Pi(r.ind, jy, jz) = SQ(Pi(r.ind, mesh->yend, jz)) / Pi(r.ind, mesh->yend - 1, jz); + } + + // Zero-gradient density + BoutReal nesheath = 0.5 * (Ne(r.ind, mesh->yend, jz) + + Ne(r.ind, mesh->yend + 1, jz)); + BoutReal tesheath = 0.5 * (Te(r.ind, mesh->yend, jz) + + Pe(r.ind, mesh->yend+1, jz) / Ne(r.ind, mesh->yend+1, jz)); + BoutReal tisheath = 0.5 * (Ti(r.ind, mesh->yend, jz) + + Pi(r.ind, mesh->yend+1, jz) / Ne(r.ind, mesh->yend+1, jz)); + + // Zero-gradient potential + BoutReal phisheath = phi(r.ind, mesh->yend, jz); + + // Ion velocity goes to the sound speed + BoutReal visheath = sqrt(tesheath + tisheath); // Sound speed outwards + + if (Vi(r.ind, mesh->yend, jz) > visheath) { + // If plasma is faster, go to plasma velocity + visheath = Vi(r.ind, mesh->yend, jz); + } + + // Sheath current + BoutReal phi_te = floor(phisheath / tesheath, 0.0); + BoutReal vesheath = + sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); + // J = n*(Vi - Ve) + BoutReal jsheath = nesheath * (visheath - vesheath); + if (nesheath < 1e-10) { + vesheath = visheath; + jsheath = 0.0; + } + + // Apply boundary condition half-way between cells + for (int jy = mesh->yend + 1; jy < mesh->LocalNy; jy++) { + // Neumann conditions + phi(r.ind, jy, jz) = phisheath; + Vort(r.ind, jy, jz) = Vort(r.ind, mesh->yend, jz); + + // Here zero-gradient Te, heat flux applied later + // This is so that heat diffusion doesn't remove (or add) additional heat + Te(r.ind, jy, jz) = Te(r.ind, mesh->yend, jz); + Ti(r.ind, jy, jz) = Ti(r.ind, mesh->yend, jz); + + // Dirichlet conditions to set flows + Vi(r.ind, jy, jz) = 2. * visheath - Vi(r.ind, mesh->yend, jz); + Ve(r.ind, jy, jz) = 2. * vesheath - Ve(r.ind, mesh->yend, jz); + Jpar(r.ind, jy, jz) = 2. * jsheath - Jpar(r.ind, mesh->yend, jz); + NVi(r.ind, jy, jz) = + 2. * nesheath * visheath - NVi(r.ind, mesh->yend, jz); + } + } + } + break; + } + default: { + throw BoutException("Sheath model %d not implemented", sheath_model); + } } } @@ -2030,7 +2425,7 @@ int Hermes::rhs(BoutReal t) { } // Collisional damping (normalised) - if (resistivity || (!electromagnetic && !FiniteElMass)) { + if (currents && (resistivity || (!electromagnetic && !FiniteElMass))) { // Need to calculate nu if electrostatic and zero electron mass nu = resistivity_multiply / (1.96 * tau_e * mi_me); @@ -2097,6 +2492,8 @@ int Hermes::rhs(BoutReal t) { } } } + + nu.applyBoundary(t); } if (thermal_conduction || sinks) { @@ -2143,11 +2540,9 @@ int Hermes::rhs(BoutReal t) { kappa_epar(r.ind, mesh->yend + 1, jz) = kappa_epar(r.ind, mesh->yend, jz); kappa_ipar(r.ind, mesh->yend + 1, jz) = kappa_ipar(r.ind, mesh->yend, jz); } - } + } } - nu.applyBoundary(t); - if (ion_viscosity) { /////////////////////////////////////////////////////////// // Ion stress tensor. Split into @@ -2187,14 +2582,15 @@ int Hermes::rhs(BoutReal t) { 0.49 * (qipar / Pilim) * (2.27 * Grad_par(log(Tilim)) - Grad_par(log(Pilim))) + 0.75 * (0.2 * SQ(qipar) - 0.085 * qisq) / (Pilim * Tilim); - + // Parallel part Pi_cipar = -0.96 * Pi * tau_i * (2. * Grad_par(Vi) + Vi * Grad_par(log(coord->Bxy))); // Could also be written as: // Pi_cipar = - // 0.96*Pi*tau_i*2.*Grad_par(sqrt(coord->Bxy)*Vi)/sqrt(coord->Bxy); - + + if (mesh->firstX()) { // First cells in X subject to boundary effects. for(int i=mesh->xstart; (i <= mesh->xend) && (i < 4); i++) { @@ -2280,7 +2676,7 @@ int Hermes::rhs(BoutReal t) { if (currents) { // Parallel wave speed increased to electron sound speed // since electrostatic & electromagnetic waves are supported - ddt(Ne) -= FV::Div_par(Ne, Ve, sqrt(mi_me) * sound_speed); + ddt(Ne) -= FV::Div_par(Ne, Ve, max_speed); } else { // Parallel wave speed is ion sound speed ddt(Ne) -= FV::Div_par(Ne, Ve, sound_speed); @@ -2290,7 +2686,7 @@ int Hermes::rhs(BoutReal t) { if (j_diamag) { // Diamagnetic drift, formulated as a magnetic drift // i.e Grad-B + curvature drift - ddt(Ne) -= FV::Div_f_v(Ne, -Telim * Curlb_B, ne_bndry_flux); + ddt(Ne) -= j_diamag_scale * FV::Div_f_v(Ne, -Telim * Curlb_B, ne_bndry_flux); } if (ramp_mesh && (t < ramp_timescale)) { @@ -2311,7 +2707,7 @@ int Hermes::rhs(BoutReal t) { } // Source - if (adapt_source) { + if (adapt_source_n) { // Add source. Ensure that sink will go to zero as Ne -> 0 Field2D NeErr = averageY(DC(Ne) - NeTarget); @@ -2328,7 +2724,7 @@ int Hermes::rhs(BoutReal t) { Sn(x, y) -= source_p * NeErr(x, y); ddt(Sn)(x, y) = -source_i * NeErr(x, y); - if (Sn(x, y) < 0.0) { + if (sources_positive && Sn(x, y) < 0.0) { Sn(x, y) = 0.0; if (ddt(Sn)(x, y) < 0.0) ddt(Sn)(x, y) = 0.0; @@ -2391,27 +2787,36 @@ int Hermes::rhs(BoutReal t) { // Note: This term is central differencing so that it balances // the corresponding compression term in the pressure equation - ddt(Vort) += Div((Pi + Pe) * Curlb_B); + ddt(Vort) += j_diamag_scale * Div((Pi + Pe) * Curlb_B); } // Advection of vorticity by ExB if (boussinesq) { TRACE("Vort:boussinesq"); // Using the Boussinesq approximation + if (j_pol_pi){ - ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, vort_bndry_flux, - poloidal_flows); + ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, vort_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) / SQ(coord->Bxy); - DelpPhi_2B2.applyBoundary("free_o2"); + // 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) / SQ(coord->Bxy); + DelpPhi_2B2.applyBoundary("free_o2"); - mesh->communicate(vEdotGradPi, DelpPhi_2B2); + mesh->communicate(vEdotGradPi, DelpPhi_2B2); - ddt(Vort) -= FV::Div_a_Laplace_perp(0.5 / SQ(coord->Bxy), vEdotGradPi); + ddt(Vort) -= FV::Div_a_Laplace_perp(0.5 / SQ(coord->Bxy), vEdotGradPi); + + ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi, vort_bndry_flux, + poloidal_flows); + } else { + // use simplified polarization term from i.e. GBS + ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(Vort, phi, vort_bndry_flux, + poloidal_flows); + } } else { // When the Boussinesq approximation is not made, @@ -2463,25 +2868,15 @@ int Hermes::rhs(BoutReal t) { if (vort_dissipation) { // Adds dissipation term like in other equations - // Maximum speed either electron sound speed or Alfven speed - Field3D max_speed = Bnorm * coord->Bxy / - sqrt(SI::mu0 * AA * SI::Mp * Nnorm * Nelim) / - Cs0; // Alfven speed (normalised by Cs0) - Field3D elec_sound = sqrt(mi_me) * sound_speed; // Electron sound speed - for (auto& i : max_speed.getRegion(RGN_ALL)) { - if (elec_sound[i] > max_speed[i]) { - max_speed[i] = elec_sound[i]; - } - - // Limit to 100x reference sound speed or light speed - BoutReal lim = BOUTMIN(100., 3e8/Cs0); - if (max_speed[i] > lim) { - max_speed[i] = lim; - } - } ddt(Vort) -= FV::Div_par(Vort, 0.0, max_speed); } + + if (phi_dissipation) { + // Adds dissipation term like in other equations, but depending on gradient of potential + // Note: Doesn't seem to need faster than sound speed + ddt(Vort) -= FV::Div_par(-phi, 0.0, sound_speed); + } } /////////////////////////////////////////////////////////// @@ -2555,24 +2950,6 @@ int Hermes::rhs(BoutReal t) { if (vepsi_dissipation) { // Adds dissipation term like in other equations - // Maximum speed either electron sound speed or Alfven speed - Field3D max_speed = Bnorm * coord->Bxy / - sqrt(SI::mu0 * AA * SI::Mp * Nnorm * Nelim) / - Cs0; // Alfven speed (normalised by Cs0) - Field3D elec_sound = sqrt(mi_me) * sound_speed; // Electron sound speed - for (auto& i : max_speed.getRegion(RGN_ALL)) { - // Maximum of Alfven or thermal electron speed - if (elec_sound[i] > max_speed[i]) { - max_speed[i] = elec_sound[i]; - } - - // Limit to 100x reference sound speed or light speed - BoutReal lim = BOUTMIN(100., 3e8/Cs0); - if (max_speed[i] > lim) { - max_speed[i] = lim; - } - } - ddt(VePsi) -= FV::Div_par(Ve - Vi, 0.0, max_speed); } } @@ -2592,8 +2969,9 @@ int Hermes::rhs(BoutReal t) { if (j_diamag) { // Magnetic drift - ddt(NVi) -= FV::Div_f_v(NVi, Tilim * Curlb_B, - ne_bndry_flux); // Grad-B, curvature drift + ddt(NVi) -= j_diamag_scale + * FV::Div_f_v(NVi, Tilim * Curlb_B, + ne_bndry_flux); // Grad-B, curvature drift } ddt(NVi) -= FV::Div_par_fvv(Ne, Vi, sound_speed, false); //FV::Div_par(NVi, Vi, sound_speed, false); @@ -2602,14 +2980,19 @@ int Hermes::rhs(BoutReal t) { if (pe_par) { ddt(NVi) -= Grad_parP(Pe + Pi); } - - if (ion_viscosity) { - TRACE("NVi:ion viscosity"); - // Poloidal flow damping + + if (ion_viscosity_par) { + TRACE("NVi:ion viscosity parallel"); + // Poloidal flow damping parallel part // The parallel part is solved as a diffusion term ddt(NVi) += 1.28 * sqrtB * FV::Div_par_K_Grad_par(Pi * tau_i / (coord->Bxy), sqrtB * Vi); + } + + if (ion_viscosity) { + TRACE("NVi:ion viscosity"); + // Poloidal flow damping if (currents) { // Perpendicular part. B32 = B^{3/2} @@ -2682,7 +3065,7 @@ int Hermes::rhs(BoutReal t) { // Parallel flow if (currents) { // Like Ne term, parallel wave speed increased - ddt(Pe) -= FV::Div_par(Pe, Ve, sqrt(mi_me) * sound_speed); + ddt(Pe) -= FV::Div_par(Pe, Ve, max_speed); } else { ddt(Pe) -= FV::Div_par(Pe, Ve, sound_speed); } @@ -2690,11 +3073,12 @@ int Hermes::rhs(BoutReal t) { if (j_diamag) { // Diamagnetic flow // Magnetic drift (curvature) divergence. - ddt(Pe) -= (5. / 3) * FV::Div_f_v(Pe, -Telim * Curlb_B, pe_bndry_flux); + ddt(Pe) -= + j_diamag_scale * (5. / 3) * FV::Div_f_v(Pe, -Telim * Curlb_B, pe_bndry_flux); // This term energetically balances diamagnetic term // in the vorticity equation - ddt(Pe) -= (2. / 3) * Pe * (Curlb_B * Grad(phi)); + ddt(Pe) -= j_diamag_scale * (2. / 3) * floor(Pe, 0.0) * (Curlb_B * Grad(phi)); } // Parallel heat conduction @@ -2733,7 +3117,8 @@ int Hermes::rhs(BoutReal t) { switch (sheath_model) { case 0: case 2: - case 3: { + case 3: + case 4: { for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Temperature and density at the sheath entrance @@ -2745,7 +3130,7 @@ int Hermes::rhs(BoutReal t) { 0.0); BoutReal nesheath = floor( 0.5 * (Ne_FA(r.ind, mesh->yend, jz) + Ne_FA(r.ind, mesh->yend + 1, jz)), - 0.0); + nesheath_floor); // Sound speed (normalised units) BoutReal Cs = sqrt(tesheath + tisheath); @@ -2754,20 +3139,23 @@ int Hermes::rhs(BoutReal t) { BoutReal q = (sheath_gamma_e - 1.5) * tesheath * nesheath * Cs; // Multiply by cell area to get power - BoutReal flux = q * (coord->J(r.ind, mesh->yend) + - coord->J(r.ind, mesh->yend + 1)) / - (sqrt(coord->g_22(r.ind, mesh->yend)) + - sqrt(coord->g_22(r.ind, mesh->yend + 1))); + BoutReal flux = q * (coord->J(r.ind, mesh->yend, jz) + + coord->J(r.ind, mesh->yend + 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->yend, jz)) + + sqrt(coord->g_22(r.ind, mesh->yend + 1, jz))); // Divide by volume of cell, and 2/3 to get pressure BoutReal power = - flux / (coord->dy(r.ind, mesh->yend) * coord->J(r.ind, mesh->yend)); + flux / (coord->dy(r.ind, mesh->yend, jz) * coord->J(r.ind, mesh->yend, jz)); sheath_dpe(r.ind, mesh->yend, jz) = -(2. / 3) * power; - wall_power(r.ind, mesh->yend) += power; + wall_power(r.ind, mesh->yend, jz) += power; } } break; } + default: { + throw BoutException("sheath_model %d not implemented", sheath_model); + } } ddt(Pe) += fromFieldAligned(sheath_dpe); } @@ -2779,7 +3167,8 @@ int Hermes::rhs(BoutReal t) { switch (sheath_model) { case 0: case 2: - case 3: { + case 3: + case 4: { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Temperature and density at the sheath entrance @@ -2791,7 +3180,7 @@ int Hermes::rhs(BoutReal t) { 0.0); BoutReal nesheath = floor(0.5 * (Ne_FA(r.ind, mesh->ystart, jz) + Ne_FA(r.ind, mesh->ystart - 1, jz)), - 0.0); + nesheath_floor); // Sound speed (normalised units) BoutReal Cs = sqrt(tesheath + tisheath); @@ -2801,20 +3190,23 @@ int Hermes::rhs(BoutReal t) { (sheath_gamma_e - 1.5) * tesheath * nesheath * Cs; // NB: positive // Multiply by cell area to get power - BoutReal flux = q * (coord->J(r.ind, mesh->ystart) + - coord->J(r.ind, mesh->ystart - 1)) / - (sqrt(coord->g_22(r.ind, mesh->ystart)) + - sqrt(coord->g_22(r.ind, mesh->ystart - 1))); + BoutReal flux = q * (coord->J(r.ind, mesh->ystart, jz) + + coord->J(r.ind, mesh->ystart - 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->ystart, jz)) + + sqrt(coord->g_22(r.ind, mesh->ystart - 1, jz))); // Divide by volume of cell, and 2/3 to get pressure - BoutReal power = flux / (coord->dy(r.ind, mesh->ystart) * - coord->J(r.ind, mesh->ystart)); + BoutReal power = flux / (coord->dy(r.ind, mesh->ystart, jz) * + coord->J(r.ind, mesh->ystart, jz)); sheath_dpe(r.ind, mesh->ystart, jz) = -(2. / 3) * power; - wall_power(r.ind, mesh->ystart) += power; + wall_power(r.ind, mesh->ystart, jz) += power; } } break; } + default: { + throw BoutException("sheath_model %d not implemented", sheath_model); + } } ddt(Pe) += fromFieldAligned(sheath_dpe); @@ -2827,7 +3219,7 @@ int Hermes::rhs(BoutReal t) { if (pe_par_p_term) { // This term balances energetically the pressure term // in Ohm's law - ddt(Pe) -= (2. / 3) * Pelim * Div_par(Ve); + ddt(Pe) -= (2. / 3) * floor(Pe, 0.0) * Div_par(Ve); } if (ramp_mesh && (t < ramp_timescale)) { ddt(Pe) += PeTarget / ramp_timescale; @@ -2860,7 +3252,7 @@ int Hermes::rhs(BoutReal t) { ////////////////////// // Sources - if (adapt_source) { + if (adapt_source_p) { // Add source. Ensure that sink will go to zero as Pe -> 0 Field2D PeErr = averageY(DC(Pe) - PeTarget); @@ -2876,7 +3268,7 @@ int Hermes::rhs(BoutReal t) { Spe(x, y) -= source_p * PeErr(x, y); ddt(Spe)(x, y) = -source_i * PeErr(x, y); - if (Spe(x, y) < 0.0) { + if (sources_positive && Spe(x, y) < 0.0) { Spe(x, y) = 0.0; if (ddt(Spe)(x, y) < 0.0) ddt(Spe)(x, y) = 0.0; @@ -2946,14 +3338,14 @@ int Hermes::rhs(BoutReal t) { if (j_diamag) { // Diamagnetic flow // Magnetic drift (curvature) divergence - ddt(Pi) -= (5. / 3) * FV::Div_f_v(Pi, Tilim * Curlb_B, pe_bndry_flux); + ddt(Pi) -= j_diamag_scale * (5. / 3) * FV::Div_f_v(Pi, Tilim * Curlb_B, pe_bndry_flux); // Compression of ExB flow // These terms energetically balances diamagnetic term // in the vorticity equation - ddt(Pi) -= (2. / 3) * Pi * (Curlb_B * Grad(phi)); + ddt(Pi) -= j_diamag_scale * (2. / 3) * Pi * (Curlb_B * Grad(phi)); - ddt(Pi) += Pi * Div((Pe + Pi) * Curlb_B); + ddt(Pi) += j_diamag_scale * Pi * Div((Pe + Pi) * Curlb_B); } if (j_par) { @@ -3054,7 +3446,8 @@ int Hermes::rhs(BoutReal t) { switch (sheath_model) { case 0: case 2: - case 3: { + case 3: + case 4: { for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Temperature and density at the sheath entrance @@ -3066,7 +3459,7 @@ int Hermes::rhs(BoutReal t) { 0.0); BoutReal nesheath = floor( 0.5 * (Ne_FA(r.ind, mesh->yend, jz) + Ne_FA(r.ind, mesh->yend + 1, jz)), - 0.0); + nesheath_floor); // Sound speed (normalised units) BoutReal Cs = sqrt(tesheath + tisheath); @@ -3075,20 +3468,23 @@ int Hermes::rhs(BoutReal t) { BoutReal q = (sheath_gamma_i - 1.5) * tisheath * nesheath * Cs; // Multiply by cell area to get power - BoutReal flux = q * (coord->J(r.ind, mesh->yend) + - coord->J(r.ind, mesh->yend + 1)) / - (sqrt(coord->g_22(r.ind, mesh->yend)) + - sqrt(coord->g_22(r.ind, mesh->yend + 1))); + BoutReal flux = q * (coord->J(r.ind, mesh->yend, jz) + + coord->J(r.ind, mesh->yend + 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->yend, jz)) + + sqrt(coord->g_22(r.ind, mesh->yend + 1, jz))); // Divide by volume of cell, and 2/3 to get pressure BoutReal power = - flux / (coord->dy(r.ind, mesh->yend) * coord->J(r.ind, mesh->yend)); + flux / (coord->dy(r.ind, mesh->yend, jz) * coord->J(r.ind, mesh->yend, jz)); sheath_dpi(r.ind, mesh->yend, jz) = -(2. / 3) * power; - wall_power(r.ind, mesh->yend) += power; + wall_power(r.ind, mesh->yend, jz) += power; } } break; } + default: { + throw BoutException("sheath_model %d not implemented", sheath_model); + } } ddt(Pi) += fromFieldAligned(sheath_dpi); } @@ -3100,7 +3496,8 @@ int Hermes::rhs(BoutReal t) { switch (sheath_model) { case 0: case 2: - case 3: { + case 3: + case 4: { for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Temperature and density at the sheath entrance @@ -3112,7 +3509,7 @@ int Hermes::rhs(BoutReal t) { 0.0); BoutReal nesheath = floor(0.5 * (Ne_FA(r.ind, mesh->ystart, jz) + Ne_FA(r.ind, mesh->ystart - 1, jz)), - 0.0); + nesheath_floor); // Sound speed (normalised units) BoutReal Cs = sqrt(tesheath + tisheath); @@ -3121,20 +3518,23 @@ int Hermes::rhs(BoutReal t) { BoutReal q = (sheath_gamma_i - 1.5) * tisheath * nesheath * Cs; // Multiply by cell area to get power - BoutReal flux = q * (coord->J(r.ind, mesh->ystart) + - coord->J(r.ind, mesh->ystart - 1)) / - (sqrt(coord->g_22(r.ind, mesh->ystart)) + - sqrt(coord->g_22(r.ind, mesh->ystart - 1))); + BoutReal flux = q * (coord->J(r.ind, mesh->ystart, jz) + + coord->J(r.ind, mesh->ystart - 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->ystart, jz)) + + sqrt(coord->g_22(r.ind, mesh->ystart - 1, jz))); // Divide by volume of cell, and 2/3 to get pressure - BoutReal power = flux / (coord->dy(r.ind, mesh->ystart) * - coord->J(r.ind, mesh->ystart)); + BoutReal power = flux / (coord->dy(r.ind, mesh->ystart, jz) * + coord->J(r.ind, mesh->ystart, jz)); sheath_dpi(r.ind, mesh->ystart, jz) = -(2. / 3) * power; - wall_power(r.ind, mesh->ystart) += power; + wall_power(r.ind, mesh->ystart, jz) += power; } } break; } + default: { + throw BoutException("sheath_model %d not implemented", sheath_model); + } } ddt(Pi) += fromFieldAligned(sheath_dpi); @@ -3143,7 +3543,7 @@ int Hermes::rhs(BoutReal t) { ////////////////////// // Sources - if (adapt_source) { + if (adapt_source_p) { // Add source. Ensure that sink will go to zero as Pe -> 0 Field2D PiErr = averageY(DC(Pi) - PiTarget); @@ -3159,7 +3559,7 @@ int Hermes::rhs(BoutReal t) { Spi(x, y) -= source_p * PiErr(x, y); ddt(Spi)(x, y) = -source_i * PiErr(x, y); - if (Spi(x, y) < 0.0) { + if (sources_positive && Spi(x, y) < 0.0) { Spi(x, y) = 0.0; if (ddt(Spi)(x, y) < 0.0) ddt(Spi)(x, y) = 0.0; @@ -3216,11 +3616,31 @@ int Hermes::rhs(BoutReal t) { if (radial_buffers) { /// Radial buffer regions - // Calculate flux sZ averages - Field2D PeDC = averageY(DC(Pe)); - Field2D PiDC = averageY(DC(Pi)); - Field2D NeDC = averageY(DC(Ne)); - Field2D VortDC = averageY(DC(Vort)); + // Calculate flux Z averages. + // This is used for both inner and outer boundaries + Field2D PeDC = DC(Pe); + Field2D PiDC = DC(Pi); + Field2D NeDC = DC(Ne); + Field2D VortDC = DC(Vort); + Field2D NViDC = DC(NVi); + Field2D VePsiDC = DC(VePsi); + + // Flux surface averages. + // In the core region it can be desirable to damp towards a flux surface average + + // First the plasma density and pressures, which should be approximately + // constant on core flux surfaces + Field2D PeInner = radial_inner_averagey ? averageY(PeDC) : PeDC; + Field2D PiInner = radial_inner_averagey ? averageY(PiDC) : PiDC; + Field2D NeInner = radial_inner_averagey ? averageY(NeDC) : NeDC; + + // Vorticity. Probably usually don't usually want to average this in Y + Field2D VortInner = radial_inner_averagey_vort ? averageY(VortDC) : VortDC; + + // Parallel flow can be damped to zero, a constant value, or allowed to vary in Y + Field2D NViInner = radial_inner_zero_nvi + ? 0.0 + : (radial_inner_averagey_nvi ? averageY(NViDC) : NViDC); if ((mesh->getGlobalXIndex(mesh->xstart) - mesh->xstart) < radial_inner_width) { // This processor contains points inside the inner radial boundary @@ -3246,23 +3666,24 @@ int Hermes::rhs(BoutReal t) { BoutReal D = radial_buffer_D * (1. - pos / radial_inner_width); for (int j = mesh->ystart; j <= mesh->yend; ++j) { - BoutReal dx = coord->dx(i, j); - BoutReal dx_xp = coord->dx(i + 1, j); - BoutReal J = coord->J(i, j); - BoutReal J_xp = coord->J(i + 1, j); - - // Calculate metric factors for radial fluxes - BoutReal rad_flux_factor = 0.25 * (J + J_xp) * (dx + dx_xp); - BoutReal x_factor = rad_flux_factor / (J * dx); - BoutReal xp_factor = rad_flux_factor / (J_xp * dx_xp); - for (int k = 0; k < ncz; ++k) { - // Relax towards constant value on flux surface - ddt(Pe)(i, j, k) -= D * (Pe(i, j, k) - PeDC(i, j)); - ddt(Pi)(i, j, k) -= D * (Pi(i, j, k) - PiDC(i, j)); - ddt(Ne)(i, j, k) -= D * (Ne(i, j, k) - NeDC(i, j)); - ddt(Vort)(i, j, k) -= D * (Vort(i, j, k) - VortDC(i, j)); - ddt(NVi)(i, j, k) -= D * NVi(i, j, k); + BoutReal dx = coord->dx(i, j, k); + BoutReal dx_xp = coord->dx(i + 1, j, k); + BoutReal J = coord->J(i, j, k); + BoutReal J_xp = coord->J(i + 1, j, k); + + // Calculate metric factors for radial fluxes + BoutReal rad_flux_factor = 0.25 * (J + J_xp) * (dx + dx_xp); + BoutReal x_factor = rad_flux_factor / (J * dx); + BoutReal xp_factor = rad_flux_factor / (J_xp * dx_xp); + + // Relax towards constant value on flux surface + ddt(Pe)(i, j, k) -= D * (Pe(i, j, k) - PeInner(i, j)); + ddt(Pi)(i, j, k) -= D * (Pi(i, j, k) - PiInner(i, j)); + ddt(Ne)(i, j, k) -= D * (Ne(i, j, k) - NeInner(i, j)); + ddt(Vort)(i, j, k) -= D * (Vort(i, j, k) - VortInner(i, j)); + ddt(NVi)(i, j, k) -= D * (NVi(i, j, k) - NViInner(i, j)); + ddt(VePsi)(i, j, k) -= D * (VePsi(i, j, k) - VePsiDC(i, j)); // Radial fluxes BoutReal f = D * (Ne(i + 1, j, k) - Ne(i, j, k)); @@ -3308,23 +3729,26 @@ int Hermes::rhs(BoutReal t) { BoutReal D = radial_buffer_D * (1. - pos / radial_outer_width); for (int j = mesh->ystart; j <= mesh->yend; ++j) { - BoutReal dx = coord->dx(i, j); - BoutReal dx_xp = coord->dx(i + 1, j); - BoutReal J = coord->J(i, j); - BoutReal J_xp = coord->J(i + 1, j); - - // Calculate metric factors for radial fluxes - BoutReal rad_flux_factor = 0.25 * (J + J_xp) * (dx + dx_xp); - BoutReal x_factor = rad_flux_factor / (J * dx); - BoutReal xp_factor = rad_flux_factor / (J_xp * dx_xp); - for (int k = 0; k < ncz; ++k) { + BoutReal dx = coord->dx(i, j, k); + BoutReal dx_xp = coord->dx(i + 1, j, k); + BoutReal J = coord->J(i, j, k); + BoutReal J_xp = coord->J(i + 1, j, k); + + // Calculate metric factors for radial fluxes + BoutReal rad_flux_factor = 0.25 * (J + J_xp) * (dx + dx_xp); + BoutReal x_factor = rad_flux_factor / (J * dx); + BoutReal xp_factor = rad_flux_factor / (J_xp * dx_xp); + ddt(Pe)(i, j, k) -= D * (Pe(i, j, k) - PeDC(i, j)); ddt(Pi)(i, j, k) -= D * (Pi(i, j, k) - PiDC(i, j)); ddt(Ne)(i, j, k) -= D * (Ne(i, j, k) - NeDC(i, j)); ddt(Vort)(i, j, k) -= D * (Vort(i, j, k) - VortDC(i, j)); - // ddt(Vort)(i,j,k) -= D*Vort(i,j,k); + ddt(NVi)(i, j, k) -= D * (NVi(i, j, k) - NViDC(i, j)); + ddt(VePsi)(i, j, k) -= D * (VePsi(i, j, k) - VePsiDC(i, j)); + // Radial fluxes + BoutReal f = D * (Vort(i + 1, j, k) - Vort(i, j, k)); ddt(Vort)(i, j, k) += f * x_factor; ddt(Vort)(i + 1, j, k) -= f * xp_factor; @@ -3381,14 +3805,14 @@ int Hermes::rhs(BoutReal t) { // Flow of neutrals inwards BoutReal flow = frecycle * flux_ion * - (coord->J(r.ind, mesh->ystart) + - coord->J(r.ind, mesh->ystart - 1)) / - (sqrt(coord->g_22(r.ind, mesh->ystart)) + - sqrt(coord->g_22(r.ind, mesh->ystart - 1))); + (coord->J(r.ind, mesh->ystart, jz) + + coord->J(r.ind, mesh->ystart - 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->ystart, jz)) + + sqrt(coord->g_22(r.ind, mesh->ystart - 1, jz))); // Rate of change of neutrals in final cell - BoutReal dndt = flow / (coord->J(r.ind, mesh->ystart) * - coord->dy(r.ind, mesh->ystart)); + BoutReal dndt = flow / (coord->J(r.ind, mesh->ystart, jz) * + coord->dy(r.ind, mesh->ystart, jz)); // Add mass, momentum and energy to the neutrals @@ -3399,7 +3823,7 @@ int Hermes::rhs(BoutReal t) { dndt * neutral_vwall * sqrt(3.5 / Tnorm)); // Power deposited onto the wall due to surface recombination - wall_power(r.ind, mesh->ystart) += (13.6 / Tnorm) * dndt; + wall_power(r.ind, mesh->ystart, jz) += (13.6 / Tnorm) * dndt; } } } @@ -3418,14 +3842,14 @@ int Hermes::rhs(BoutReal t) { 0.5 * (Ve(r.ind, mesh->yend, jz) + Ve(r.ind, mesh->yend + 1, jz)); // Flow of neutrals inwards - BoutReal flow = flux_ion * (coord->J(r.ind, mesh->yend) + - coord->J(r.ind, mesh->yend + 1)) / - (sqrt(coord->g_22(r.ind, mesh->yend)) + - sqrt(coord->g_22(r.ind, mesh->yend + 1))); + BoutReal flow = flux_ion * (coord->J(r.ind, mesh->yend, jz) + + coord->J(r.ind, mesh->yend + 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->yend, jz)) + + sqrt(coord->g_22(r.ind, mesh->yend + 1, jz))); // Rate of change of neutrals in final cell BoutReal dndt = - flow / (coord->J(r.ind, mesh->yend) * coord->dy(r.ind, mesh->yend)); + flow / (coord->J(r.ind, mesh->yend, jz) * coord->dy(r.ind, mesh->yend, jz)); // Add mass, momentum and energy to the neutrals @@ -3436,7 +3860,7 @@ int Hermes::rhs(BoutReal t) { -dndt * neutral_vwall * sqrt(3.5 / Tnorm)); // Power deposited onto the wall due to surface recombination - wall_power(r.ind, mesh->yend) += (13.6 / Tnorm) * dndt; + wall_power(r.ind, mesh->yend, jz) += (13.6 / Tnorm) * dndt; } } } @@ -3593,7 +4017,7 @@ int Hermes::precon(BoutReal t, BoutReal gamma, BoutReal delta) { static std::unique_ptr inv; if (!inv) { // Initialise parallel inversion class - inv = InvertPar::Create(); + auto inv = InvertPar::create(); inv->setCoefA(1.0); } if (thermal_conduction) { diff --git a/hermes-2.hxx b/hermes-2.hxx index ad34aa7..a43cb30 100644 --- a/hermes-2.hxx +++ b/hermes-2.hxx @@ -48,6 +48,8 @@ private: // Equilibrium current Field2D Jpar0; + BoutReal nesheath_floor; // Density floor used in sheath boundary conditions + // Evolving variables Field3D Ne; // Electron density Field3D Pe, Pi; // Electron and Ion pressures @@ -101,7 +103,10 @@ private: bool FiniteElMass; // Finite Electron Mass bool j_diamag; // Diamagnetic current: Vort <-> Pe + FieldGeneratorPtr j_diamag_scale_generator; // Time-varying diamagnetic current scaling + BoutReal j_diamag_scale; // Diamagnetic current scaling factor. bool j_par; // Parallel current: Vort <-> Psi + bool j_pol_pi; // Explicit Pi in polarization current bool parallel_flow; bool parallel_flow_p_term; // Vi advection terms in Pe, Pi bool pe_par; // Parallel pressure gradient: Pe <-> Psi @@ -110,6 +115,7 @@ private: bool thermal_force; // Force due to temperature gradients bool electron_viscosity; // Electron parallel viscosity bool ion_viscosity; // Ion viscosity + bool ion_viscosity_par; // Parallel part of ion viscosity bool electron_neutral; // Include electron-neutral collisions in resistivity bool ion_neutral; // Include ion-neutral collisions in ion collision time bool poloidal_flows; // Include y derivatives in diamagnetic and ExB drifts @@ -129,8 +135,6 @@ private: bool ion_velocity; // Include Vi terms bool phi3d; // Use a 3D solver for phi - - bool staggered; // Use staggered differencing along B bool boussinesq; // Use a fixed density (Nnorm) in the vorticity equation @@ -142,6 +146,10 @@ private: int radial_inner_width; // Number of points in the inner radial buffer int radial_outer_width; // Number of points in the outer radial buffer BoutReal radial_buffer_D; // Diffusion in buffer region + bool radial_inner_averagey; // Average Ne, Pe, Pi fields in Y in inner radial buffer + bool radial_inner_averagey_vort; // Average vorticity in Y in inner buffer + bool radial_inner_averagey_nvi; // Average NVi in Y in inner buffer + bool radial_inner_zero_nvi; // Damp NVi towards zero in inner buffer BoutReal resistivity_boundary; // Value of nu in boundary layer int resistivity_boundary_width; // Width of radial boundary @@ -163,7 +171,7 @@ private: // Fix density in SOL bool sol_fix_profiles; std::shared_ptr sol_ne, sol_te; // Generating functions - + // Output switches for additional information bool verbose; // Outputs additional fields, mainly for debugging bool output_ddt; // Output time derivatives @@ -180,6 +188,7 @@ private: BoutReal floor_num_cs; // Apply a floor to the numerical sound speed bool vepsi_dissipation; // Dissipation term in VePsi equation bool vort_dissipation; // Dissipation term in Vorticity equation + bool phi_dissipation; // Dissipation term in Vorticity equation, depending on phi // Sources and profiles @@ -187,7 +196,9 @@ private: BoutReal ramp_timescale; // Length of time for the initial ramp Field2D NeTarget, PeTarget, PiTarget; // For adaptive sources - bool adapt_source; // Use a PI controller to feedback profiles + bool adapt_source_p; // Use a PI controller to feedback pressure profiles + bool adapt_source_n; // Use a PI controller to feedback density profiles + bool sources_positive; // Ensure sources > 0 bool core_sources; // Sources only in the core bool energy_source; // Add the same amount of energy to each particle BoutReal source_p, source_i; // Proportional-Integral controller @@ -196,7 +207,7 @@ private: bool density_inflow; // Does incoming density have momentum? bool source_vary_g11; // Multiply source by g11 - Field2D g11norm; + Coordinates::FieldMetric g11norm; // Boundary fluxes @@ -219,7 +230,8 @@ private: // Electromagnetic solver for finite electron mass case bool split_n0_psi; // Split the n=0 component of Apar (psi)? //Laplacian *aparSolver; - std::unique_ptr aparSolver; + + std::unique_ptr aparSolver{nullptr}; LaplaceXY *aparXY; // Solves n=0 component Field2D psi2D; // Axisymmetric Psi @@ -228,13 +240,18 @@ private: bool split_n0; // Split solve into n=0 and n~=0? LaplaceXY *laplacexy; // Laplacian solver in X-Y (n=0) Field2D phi2D; // Axisymmetric phi + + bool phi_boundary_relax; ///< Relax the boundary towards Neumann? + BoutReal phi_boundary_timescale; ///< Relaxation timescale + BoutReal phi_boundary_last_update; ///< The last time the boundary was updated bool newXZsolver; - std::unique_ptr phiSolver; // Old Laplacian in X-Z - std::unique_ptr newSolver; // New Laplacian in X-Z + + std::unique_ptr phiSolver{nullptr}; // Old Laplacian in X-Z + std::unique_ptr newSolver{nullptr}; // New Laplacian in X-Z // Mesh quantities - Field2D B32, sqrtB; + Coordinates::FieldMetric B32, sqrtB; }; /// Fundamental constants diff --git a/input_sources.py b/input_sources.py new file mode 100755 index 0000000..d5261da --- /dev/null +++ b/input_sources.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python +# +# Run this with a directory to print the input sources +# +# $ ./input_source.py data/ + +import numpy as np + +from boutdata import collect + + +def total_sources(path, tind=-1): + """Return total sources integrated over the domain in SI units: + - Volume [m^3] + - Particles [s^-1] + - Power into electrons and ions [Watts] + + Inputs + ------ + path The path to the data files + tind Time index + """ + + # Normalisation + qe = 1.602e-19 # Electron charge [C] + Nnorm = collect("Nnorm", path=path) # Reference density [m^-3] + Tnorm = collect("Tnorm", path=path) # Reference temperature [eV] + rho_s0 = collect("rho_s0", path=path) # Reference length [m] + wci = collect("Omega_ci", path=path) # Reference frequency [s^-1] + + # Metrics + J = collect("J", path=path) + dx = collect("dx", path=path) + dy = collect("dy", path=path) + dz = collect("dz", path=path) + + dV = J * dx * dy * dz # Normalised cell volume + + dV_si = dV * rho_s0 ** 3 # Cell volume [m^3] as a function of [x,y] + + result = {} + + # Read and convert units for each source + # Note the 3/2 factor multiplying the pressure source to convert to internal energy + for source_name, norm in [ + ("NeSource", Nnorm * wci), + ("PeSource", (3.0 / 2) * qe * Nnorm * Tnorm * wci), + ("PiSource", (3.0 / 2) * qe * Nnorm * Tnorm * wci), + ]: + # Read source, select single time, exclude X guard cells, convert to SI units + source = norm * collect(source_name, path=path, tind=tind)[-1, 2:-2, :, :] + + # Get number of Z points + nz = source.shape[-1] + + # Multiply by cell volume [m^3] + for z in range(nz): + source[:, :, z] *= dV_si[2:-2, :] + + # Sum over all cells + result[source_name] = np.sum(source) + + result["volume"] = np.sum(dV_si[2:-2, :]) * nz # Exclude guard cells + return result + + +if __name__ == "__main__": + + import sys + + if len(sys.argv) != 2: + # Print usage information + print("Usage: {0} path\n e.g. {0} data".format(sys.argv[0])) + sys.exit(1) + + path = sys.argv[1] + sources = total_sources(path) + + print("Volume of the simulation: {} m^3".format(sources["volume"])) + print("Particle source: {} /s".format(sources["NeSource"])) + print("Electron power: {} W".format(sources["PeSource"])) + print("Ion power: {} W".format(sources["PiSource"])) diff --git a/linear-device/README.md b/linear-device/README.md new file mode 100644 index 0000000..6227459 --- /dev/null +++ b/linear-device/README.md @@ -0,0 +1,19 @@ +Linear device simulations +========================= + +Simulations without magnetic curvature. The diamagnetic current +term is usually turned off, because it has zero divergence. + +Annulus +------- + +Sets up a simulation where the Z coordinate is the azimuthal angle, +and X is a radial coordinate. This coordinate system has a pole +at the axis, and so X covers a range of radius with a hole in the middle. + +The parameters are set up to be approximately LAPD-like. + +To run the simulation: + + $ mpirun -np 16 ./hermes-2 -d annulus + diff --git a/linear-device/annulus/BOUT.inp b/linear-device/annulus/BOUT.inp new file mode 100644 index 0000000..4cf903a --- /dev/null +++ b/linear-device/annulus/BOUT.inp @@ -0,0 +1,301 @@ +# +# Some values taken from +# https://doi.org/10.1063/1.4759010 +# https://aip.scitation.org/doi/10.1063/1.3527987 +# + +NOUT = 1000 # number of output time-steps +TIMESTEP = 50 # time between outputs + +[restart] +init_missing = true # Set missing variables to 0 if missing in restarts + +[mesh] + +nx = 64 # Radial resolution including 4 guard cells +ny = 16 # Parallel direction +nz = 64 # number of points in azimuthal direction + +length = 17 # length of machine in meters +Rmin = 0.1 # minimum radius in meters +Rmax = 0.4 # maximum radius in meters + +Bxy = 0.1 # Magnetic field strength [T] + +# The following choices make a Clebsch coordinate system where +# x is a radial flux coordinate +# y is a parallel angle coordinate (0 -> 2π) +# z is azimuthal angle (0 -> 2π) +# +# Note: In input expressions, +# x is normalised from 0->1, y and z from 0->2π, + +Bpxy = Bxy +Btxy = 0 +hthe = length / (2π) +Rxy = Rmin + (Rmax - Rmin) * x # Radius from axis. Note: Here x is from 0->1 +sinty = 0 # Integrated shear + +dr = (Rmax - Rmin) / (nx - 4) +dx = Bpxy * Rxy * dr # Radial flux coordinate +dy = 2π / ny # Angle 0 -> 2π +dz = 2π / nz # Azimuthal angle + +ixseps1 = -1 # This line and the one below will add y boundaries +ixseps2 = -1 # + +paralleltransform = identity + +Ne0 = 1e-3 * exp(-x^2) # Starting density profile [x10^20 m^-3] +Te0 = 1 # Starting electron temperature [eV] +Ti0 = 0.1 # Starting ion temperate [eV] + +################################################## +# Derivative methods + +[mesh:ddx] + +first = C2 +second = C2 +upwind = W3 + +[mesh:ddy] + +first = C2 +second = C2 +upwind = W3 + +[mesh:ddz] + +first = FFT +second = FFT +upwind = W3 + +################################################### +# Time-integration solver + +[solver] +# Note: If evolving neutrals, need preconditioning +type = cvode +use_precon = true + +ATOL = 1.0e-10 # absolute tolerance +RTOL = 1.0e-5 # relative tolerance +mxstep = 1000000 # Maximum internal steps per output + +cvode_max_order = 2 +cvode_stability_limit_detection = true + +################################################## +# Electrostatic potential solver + +[phiSolver] +inner_boundary_flags = 1 # Zero gradient +outer_boundary_flags = 16 # INVERT_SET, setting outer boundary + +[laplacexy] # 2D solver in X-Y +pctype = hypre # Preconditioner + +atol = 1e-12 +rtol = 1e-8 + +# Note: These boundary flags should be consistent with phiSolver +core_bndry_dirichlet = false +pf_bndry_dirichlet = true +y_bndry_dirichlet = false + +include_y_derivs = true # Include the y components in the inversion + +[aparSolver] +inner_boundary_flags = 0 +outer_boundary_flags = 0 + +################################################### +# Model settings + +[Hermes] + +ramp_mesh = false + +####################### +# Output variables +output_ddt = false # Output time derivatives +verbose = true # Output more derived quantities + +####################### +# Numerical dissipation + +vepsi_dissipation = true # Parallel dissipation on Ve-Vi +vort_dissipation = false # Parallel dissipation of vorticity oscillations +phi_dissipation = true # Parallel dissipation of potential oscillations + +radial_buffers = true # Turn on dissipation close to boundaries +radial_inner_width = 4 # Number of cells on inner boundary +radial_outer_width = 4 # Number of cell on outer boundary. < 0 -> no buffer +radial_buffer_D = 1 # Damping coefficient in boundary + +####################### +# Flux limiters +kappa_limit_alpha = 0.2 # SOLPS style heat flux limiter +eta_limit_alpha = 0.5 # SOLPS style viscosity limiter + +####################### +# Electric field and Ohm's law +electromagnetic = true # Electromagnetic? Otherwise electrostatic +FiniteElMass = true # Finite electron mass? (false -> Zero electron mass) + +# Electrostatic potential +split_n0 = false # Solve n=0 electrostatic potential separately? +split_n0_psi = false # Solve n=0 Apar component separately? + +phi_boundary_relax = true # Relax the boundary towards zero gradient +phi_boundary_timescale = 1e-5 # Timescale for relaxation [seconds] + +####################### +# Terms included in the model + +# NOTE: all currents switched off for fluid run +j_diamag = false # Diamagnetic current: Vort <-> Pe +j_par = true # Parallel current: Vort <-> Psi + +pe_par = true # Parallel pressure gradient: Pe <-> Psi +resistivity = true # Resistivity: Psi -> Pe +thermal_flux = true +thermal_force = true # Parallel thermal force in Ohm's law +electron_viscosity = true # Electron viscosity in Ohm's law +ion_viscosity = false # Viscous terms in vorticity equation +ion_viscosity_par = true # Parallel diffusion of ion momentum +thermal_conduction = true # Parallel heat conduction +poloidal_flows = false # Include ExB flows in the X-Y plane? +ion_velocity = true # Evolve ion parallel velocity? +boussinesq = true # Use Boussinesq approximation? +ion_neutral = false # Include ion-neutral collisions in tau_i? +neutral_friction = false # Friction between plasma and neutrals in vorticity equation + +frecycle = 0.9 # Neutral recycling fraction at targets + +carbon_fraction = -1 # Fixed fraction carbon impurity. < 0 means none + +####################### +# Transport coefficients +classical_diffusion = false # Include collisional diffusion? + +anomalous_D = -1 # Anomalous density diffusion [m^2/s] +anomalous_chi = -1 # Anomalous thermal diffusion [m^2/s] +anomalous_nu = -1 # Anomalous viscosity [m^2/s] + +####################### +# Radial boundary fluxes +ne_bndry_flux = true # Allow radial flows of particles through boundary +pe_bndry_flux = true # Allow radial flows of energy through boundary +vort_bndry_flux = false # Allow radial flows of vorticity through boundary + +####################### +# Plasma sheath +sheath_model = 2 # 0 = Bohn, 1 = Loizu, 2 = Bohm + free density +sheath_yup = true +sheath_ydown = true +sheath_gamma_e = 4 # Electron sheath heat transmission +sheath_gamma_i = 2.5 # Ion sheath heat transmission + +####################### +# Sources +core_sources = false # Only sources in the core +adapt_source = false # Feedback on profiles (PI controller) +energy_source = true # Constant power per particle (at fixed x). False -> Constant power per volume +source_p = 1e-2 # Rate based on current error (1/timescale) +source_i = 1e-6 # Integral of error (1/timescale^2) +source_vary_g11 = false # Vary source in poloidal angle to better match radial transport + +####################### +# Normalisation factors +Nnorm = 1e18 +Tnorm = 1 +Bnorm = 0.1 +AA = 4 # Atomic mass. 1 = Hydrogen, 2 = Deuterium + +####################### +# Neutrals +[neutral] +type = none # Neutral model: none, diffusion2d, recycling, fullvelocity, mixed +viscosity = 1 # Dynamic viscosity +bulk = 0 # Bulk (volume) viscosity +conduction = 1 +neutral_gamma = 0.0 # Sheath heat transmission + +nn_floor = 1e-2 # Floor applied when calculating Vn = NVn / Nn +low_n_equilibriate = -1e-4 # If n < nn_floor, evolve Tn and Vn towards plasma values + +[All] +scale = 0.0 + +bndry_all = neumann_o2 +bndry_xin = neumann_o2 +bndry_xout = neumann_o2 + +[Ne] # Electron density +scale = 1 +function = 1e-3*(mixmode(z) + mixmode(4*z - x)) +source = 6.7e2*exp(-(x/0.3)^2) # Particle source + +[Vort] +function = 0 + +[VePsi] # Ve + 0.5*beta_e*mi_me*psi +bndry_xin = dirichlet_o2 +bndry_xout = dirichlet_o2 + +[Pe] # Electron pressure +function = 1e-3 +source = 8e3*exp(-(x/0.25)^2) # Power into electrons + +[Pi] # Ion pressure +function = 1e-3 +source = 0 # Power into ions + +[Ve] + +[phi] +# Radial boundaries determined by Laplacian inversion +bndry_all = none +bndry_yup = free_o3 +bndry_ydown = neumann + +[Nn] # Neutral density + +scale = 1e-2 +function = 1 +bndry_all = neumann_o2 + +[Vn] # Neutral parallel velocity + +scale = 0 +function = 0 + +bndry_all = dirichlet_o2 + +[Vn_x] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_y] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_z] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Pn] # Neutral pressure + +scale = 1e-5 +function = 1 + +bndry_all = neumann_o2 + +[NVn] +bndry_all = dirichlet_o2 + diff --git a/linear-device/hermes-2 b/linear-device/hermes-2 new file mode 120000 index 0000000..d88e98a --- /dev/null +++ b/linear-device/hermes-2 @@ -0,0 +1 @@ +../hermes-2 \ No newline at end of file diff --git a/loadmetric.cxx b/loadmetric.cxx index cf7f7de..77bd402 100644 --- a/loadmetric.cxx +++ b/loadmetric.cxx @@ -35,7 +35,7 @@ void LoadMetric(BoutReal Lnorm, BoutReal Bnorm) { coord->Bxy /= Bnorm; // Calculate metric components - if (Options::root()["mesh"]["paralleltransform"].as() == "shifted") { + if (Options::root()["mesh"]["paralleltransform:type"].as() == "shifted") { sinty = 0.0; // I disappears from metric } diff --git a/mixed.cxx b/mixed.cxx index 8d841da..d0cd8d2 100644 --- a/mixed.cxx +++ b/mixed.cxx @@ -1,9 +1,9 @@ #include "mixed.hxx" +#include #include #include -#include #include "div_ops.hxx" @@ -15,26 +15,49 @@ NeutralMixed::NeutralMixed(Solver *solver, Mesh *UNUSED(mesh), Options &options) solver->add(Pn, "Pn"); solver->add(NVn, "NVn"); - OPTION(options, sheath_ydown, true); - OPTION(options, sheath_yup, true); - + sheath_ydown = options["sheath_ydown"] + .doc("Enable wall boundary conditions at ydown") + .withDefault(true); + + sheath_yup = options["sheath_yup"] + .doc("Enable wall boundary conditions at yup") + .withDefault(true); + neutral_gamma = options["neutral_gamma"] .doc("Heat flux to the wall q = neutral_gamma * n * T * cs") .withDefault(5. / 4); nn_floor = options["nn_floor"] - .doc("A minimum density used when dividing NVn by Nn. Normalised units.") - .withDefault(1e-4); - + .doc("A minimum density used when dividing NVn by Nn. " + "Normalised units.") + .withDefault(1e-4); + + precondition = options["precondition"] + .doc("Enable preconditioning in neutral model?") + .withDefault(true); + + if (precondition) { + inv = std::unique_ptr(Laplacian::create(&options["precon_laplace"])); + + inv->setInnerBoundaryFlags(INVERT_DC_GRAD | INVERT_AC_GRAD); + inv->setOuterBoundaryFlags(INVERT_DC_GRAD | INVERT_AC_GRAD); + + inv->setCoefA(1.0); + } + // Optionally output time derivatives - if (options["output_ddt"].doc("Save derivatives to output?").withDefault(false)) { + if (options["output_ddt"] + .doc("Save derivatives to output?") + .withDefault(false)) { SAVE_REPEAT(ddt(Nn), ddt(Pn), ddt(NVn)); } - if (options["diagnose"].doc("Save additional diagnostic fields?").withDefault(false)) { + if (options["diagnose"] + .doc("Save additional diagnostic fields (Dnn) ?") + .withDefault(false)) { SAVE_REPEAT(Dnn); } - + Dnn = 0.0; // Neutral gas diffusion S = 0; @@ -54,7 +77,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, NVn.clearParallelSlices(); Coordinates *coord = mesh->getCoordinates(); - + ////////////////////////////////////////////////////// // 3D model, diffusive in X-Z, fluid in Y // @@ -67,7 +90,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, // Nnlim Used where division by neutral density is needed Field3D Nnlim = floor(Nn, nn_floor); Field3D Tn = Pn / Nn; - //Tn = floor(Tn, 0.01 / Tnorm); + // Tn = floor(Tn, 0.01 / Tnorm); Vn = NVn / Nn; Field3D Vnlim = NVn / Nnlim; // Neutral parallel velocity @@ -101,7 +124,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, // Zero-gradient pressure Pn(r.ind, mesh->ystart - 1, jz) = Pn(r.ind, mesh->ystart, jz); - + // No flow into wall Vn(r.ind, mesh->ystart - 1, jz) = -Vn(r.ind, mesh->ystart, jz); Vnlim(r.ind, mesh->ystart - 1, jz) = -Vnlim(r.ind, mesh->ystart, jz); @@ -132,7 +155,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, // Zero-gradient pressure Pn(r.ind, mesh->yend + 1, jz) = Pn(r.ind, mesh->yend, jz); - + // No flow into wall Vn(r.ind, mesh->yend + 1, jz) = -Vn(r.ind, mesh->yend, jz); Vnlim(r.ind, mesh->yend + 1, jz) = -Vnlim(r.ind, mesh->yend, jz); @@ -152,7 +175,7 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, BoutReal neutral_lmax = 0.1 / Lnorm; Field3D Rnn = Nn * sqrt(Tn) / neutral_lmax; // Neutral-neutral collisions Dnn = Pnlim / (Riz + Rcx + Rnn); - + mesh->communicate(Dnn); Dnn.clearParallelSlices(); Dnn.applyBoundary("dirichlet_o2"); @@ -163,32 +186,32 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, Field3D DnnPn = Dnn * Pn; DnnPn.applyBoundary("dirichlet_o2"); Field3D DnnNn = Dnn * Nn; - DnnNn.applyBoundary("dirichlet_o2"); + DnnNn.applyBoundary("dirichlet_o2"); Field3D DnnNVn = Dnn * NVn; DnnNVn.applyBoundary("dirichlet_o2"); - + 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); - DnnNVn(r.ind, mesh->ystart-1, jz) = -DnnNVn(r.ind, mesh->ystart, 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); + DnnNVn(r.ind, mesh->ystart - 1, jz) = -DnnNVn(r.ind, mesh->ystart, jz); } } } - + if (sheath_yup) { 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); - DnnNVn(r.ind, mesh->yend+1, jz) = -DnnNVn(r.ind, mesh->yend, 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); + DnnNVn(r.ind, mesh->yend + 1, jz) = -DnnNVn(r.ind, mesh->yend, jz); } } } - + // Logarithms used to calculate perpendicular velocity // V_perp = -Dnn * ( Grad_perp(Nn)/Nn + Grad_perp(Tn)/Tn ) // @@ -202,15 +225,14 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, // Sound speed appearing in Lax flux for advection terms Field3D sound_speed = sqrt(Tn * (5. / 3)); - + ///////////////////////////////////////////////////// // Neutral density TRACE("Neutral density"); - ddt(Nn) = - -FV::Div_par(Nn, Vnlim, sound_speed) // Advection - + S // Source from recombining plasma - + FV::Div_a_Laplace_perp(DnnNn, logPnlim) // Perpendicular diffusion + ddt(Nn) = -FV::Div_par(Nn, Vnlim, sound_speed) // Advection + + S // Source from recombining plasma + + FV::Div_a_Laplace_perp(DnnNn, logPnlim) // Perpendicular diffusion ; ///////////////////////////////////////////////////// @@ -224,30 +246,28 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, // eta_n = (2. / 5) * kappa_n; ddt(NVn) = - -FV::Div_par(NVn, Vnlim, sound_speed) // Momentum flow - + F // Friction with plasma - - Grad_par(Pnlim) // Pressure gradient - + FV::Div_a_Laplace_perp(DnnNVn, logPnlim) // Perpendicular diffusion - + FV::Div_a_Laplace_perp((2./5)*DnnNn, Vn) // Perpendicular viscosity - + FV::Div_par_K_Grad_par((2./5)*DnnNn, Vn) // Parallel viscosity + -FV::Div_par(NVn, Vnlim, sound_speed) // Momentum flow + + F // Friction with plasma + - Grad_par(Pnlim) // Pressure gradient + + FV::Div_a_Laplace_perp(DnnNVn, logPnlim) // Perpendicular diffusion + + FV::Div_a_Laplace_perp((2. / 5) * DnnNn, Vn) // Perpendicular viscosity + + FV::Div_par_K_Grad_par((2. / 5) * DnnNn, Vn) // Parallel viscosity ; - Fperp = Rcx + Rrc; ///////////////////////////////////////////////////// // Neutral pressure TRACE("Neutral pressure"); - ddt(Pn) = - -FV::Div_par(Pn, Vnlim, sound_speed) // Advection - - (2. / 3) * Pnlim * Div_par(Vn) // Compression - + (2. / 3) * Qi // Energy exchange with ions - + FV::Div_a_Laplace_perp(DnnPn, logPnlim) // Perpendicular diffusion - + FV::Div_a_Laplace_perp(DnnNn, Tn) // Conduction - + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction + ddt(Pn) = -FV::Div_par(Pn, Vnlim, sound_speed) // Advection + - (2. / 3) * Pnlim * Div_par(Vn) // Compression + + (2. / 3) * Qi // Energy exchange with ions + + FV::Div_a_Laplace_perp(DnnPn, logPnlim) // Perpendicular diffusion + + FV::Div_a_Laplace_perp(DnnNn, Tn) // Conduction + + FV::Div_par_K_Grad_par(DnnNn, Tn) // Parallel conduction ; - + ////////////////////////////////////////////////////// // Boundary condition on fluxes TRACE("Neutral boundary fluxes"); @@ -275,15 +295,15 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, // gamma * n * T * cs BoutReal q = neutral_gamma * Nnout * Tnout * sqrt(Tnout); // Multiply by cell area to get power - BoutReal heatflux = q * (coord->J(r.ind, mesh->ystart) + - coord->J(r.ind, mesh->ystart - 1)) / - (sqrt(coord->g_22(r.ind, mesh->ystart)) + - sqrt(coord->g_22(r.ind, mesh->ystart - 1))); + BoutReal heatflux = q * (coord->J(r.ind, mesh->ystart, jz) + + coord->J(r.ind, mesh->ystart - 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->ystart, jz)) + + sqrt(coord->g_22(r.ind, mesh->ystart - 1, jz))); // Divide by volume of cell, and multiply by 2/3 to get pressure ddt(Pn)(r.ind, mesh->ystart, jz) -= - (2. / 3) * heatflux / - (coord->dy(r.ind, mesh->ystart) * coord->J(r.ind, mesh->ystart)); + (2. / 3) * heatflux / + (coord->dy(r.ind, mesh->ystart, jz) * coord->J(r.ind, mesh->ystart, jz)); } } } @@ -312,14 +332,14 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, BoutReal q = neutral_gamma * Nnout * Tnout * sqrt(Tnout); // Multiply by cell area to get power BoutReal heatflux = - q * (coord->J(r.ind, mesh->yend) + coord->J(r.ind, mesh->yend + 1)) / - (sqrt(coord->g_22(r.ind, mesh->yend)) + - sqrt(coord->g_22(r.ind, mesh->yend + 1))); + q * (coord->J(r.ind, mesh->yend, jz) + coord->J(r.ind, mesh->yend + 1, jz)) / + (sqrt(coord->g_22(r.ind, mesh->yend, jz)) + + sqrt(coord->g_22(r.ind, mesh->yend + 1, jz))); // Divide by volume of cell, and multiply by 2/3 to get pressure ddt(Pn)(r.ind, mesh->yend, jz) -= (2. / 3) * heatflux / - (coord->dy(r.ind, mesh->yend) * coord->J(r.ind, mesh->yend)); + (coord->dy(r.ind, mesh->yend, jz) * coord->J(r.ind, mesh->yend, jz)); } } } @@ -333,3 +353,19 @@ void NeutralMixed::update(const Field3D &Ne, const Field3D &Te, } } } + +void NeutralMixed::precon(BoutReal UNUSED(t), BoutReal gamma, + BoutReal UNUSED(delta)) { + if (!precondition) { + return; + } + + // Neutral gas diffusion + // Solve (1 - gamma*Dnn*Delp2)^{-1} + + inv->setCoefD(-gamma * Dnn); + + ddt(Nn) = inv->solve(ddt(Nn)); + ddt(NVn) = inv->solve(ddt(NVn)); + ddt(Pn) = inv->solve(ddt(Pn)); +} diff --git a/mixed.hxx b/mixed.hxx index 88cdb8e..68926b0 100644 --- a/mixed.hxx +++ b/mixed.hxx @@ -7,6 +7,9 @@ #define __NEUTRAL_MIXED_H__ #include "neutral-model.hxx" +#include + +#include class NeutralMixed : public NeutralModel { public: @@ -22,6 +25,8 @@ public: Field3D getDensity() override {return Nn;} + /// Preconditioner function + void precon(BoutReal t, BoutReal gamma, BoutReal delta); private: Field3D Nn, Pn, NVn; // Density, pressure and parallel momentum @@ -34,6 +39,9 @@ private: BoutReal neutral_gamma; // Heat transmission for neutrals BoutReal nn_floor; // Minimum Nn used when dividing NVn by Nn to get Vn. + + bool precondition {true}; // Enable preconditioner? + std::unique_ptr inv; // Laplacian inversion used for preconditioning }; diff --git a/neutral-model.cxx b/neutral-model.cxx index 8722b4f..94d04f7 100644 --- a/neutral-model.cxx +++ b/neutral-model.cxx @@ -105,9 +105,9 @@ void NeutralModel::neutral_rates( Nn_R = 0.0; // Jacobian (Cross-sectional area) - BoutReal J_C = coord->J(i, j), - J_L = 0.5 * (coord->J(i, j - 1) + coord->J(i, j)), - J_R = 0.5 * (coord->J(i, j) + coord->J(i, j + 1)); + BoutReal J_C = coord->J(i, j, k), + J_L = 0.5 * (coord->J(i, j - 1, k) + coord->J(i, j, k)), + J_R = 0.5 * (coord->J(i, j, k) + coord->J(i, j + 1, k)); /////////////////////////////////////////// // Charge exchange diff --git a/recycling.cxx b/recycling.cxx index 6e8f875..739a651 100644 --- a/recycling.cxx +++ b/recycling.cxx @@ -27,7 +27,11 @@ void NeutralRecycling::update(const Field3D &Ne, const Field3D &Te, TRACE("NeutralRecycling::update"); Coordinates *coord = mesh->getCoordinates(); - + + Field2D dx2D = DC(coord->dx); + Field2D dy2D = DC(coord->dy); + Field2D g_222D = DC(coord->g_22); + Field2D J2D = DC(coord->J); // Lower limit for neutral density (mainly for first time through when Nn = 0) Nn = floor(Nn, 1e-8); @@ -44,7 +48,7 @@ void NeutralRecycling::update(const Field3D &Ne, const Field3D &Te, static bool first_time = true; if (first_time) { Field2D ll; - ll = CumSumY2D(hthe * coord->dy / Lmax, true); + ll = CumSumY2D(hthe * dy2D / Lmax, true); Nn = max(Nelim) * exp(-ll); first_time = false; } @@ -76,11 +80,11 @@ void NeutralRecycling::update(const Field3D &Ne, const Field3D &Te, // number of particles lost per dt (flux out times perp volume) [t^-1] nlost = bcast_lasty(fluxout * 0.5 * - (coord->J(i, mesh->yend) * coord->dx(i, mesh->yend) * - coord->dz / sqrt(coord->g_22(i, mesh->yend)) + - coord->J(i, mesh->yend + 1) * - coord->dx(i, mesh->yend + 1) * coord->dz / - sqrt(coord->g_22(i, mesh->yend + 1)))); + (coord->J(i, mesh->yend, k) * coord->dx(i, mesh->yend, k) * + coord->dz(i, mesh->yend, k) / sqrt(coord->g_22(i, mesh->yend, k)) + + coord->J(i, mesh->yend + 1, k) * + coord->dx(i, mesh->yend + 1, k) * coord->dz(i, mesh->yend, k) / + sqrt(coord->g_22(i, mesh->yend + 1, k)))); // Integrate ionization rate over volume to get volume loss rate (simple // integration using trap rule) @@ -88,8 +92,8 @@ void NeutralRecycling::update(const Field3D &Ne, const Field3D &Te, for (int j = mesh->ystart; j <= mesh->yend; j++) { sigma_iz = hydrogen.ionisation(Te(i, j, k) * Tnorm) * Nnorm / Fnorm; // ionization rate [d^3]/[t] - BoutReal dV = coord->J(i, j) * coord->dx(i, j) * coord->dy(i, j) * - coord->dz; // volume element + BoutReal dV = coord->J(i, j, k) * coord->dx(i, j, k) * coord->dy(i, j, k) * + coord->dz (i, j, k); // volume element nnexp += Nelim(i, j, k) * sigma_iz * exp(-lambda_int(i, j, k)) * dV; // full integral of density source [d^3]/[t] } diff --git a/slab3d/BOUT.inp b/slab3d/BOUT.inp index 16a2996..1e871d9 100644 --- a/slab3d/BOUT.inp +++ b/slab3d/BOUT.inp @@ -1,16 +1,13 @@ +nout = 200 # number of output time-steps +timestep = 50 # time between outputs -NOUT = 50 -TIMESTEP = 10 +ShiftWithoutTwist = true -TwistShift = false # use twist-shift condition? -ballooning = true -shiftinitial = false -ShiftWithoutTwist=true +[input] +error_on_unused_options=false -# ZPERIOD=10 - -[restart] -init_missing=true +[restart_files] +init_missing = true # Set missing variables to 0 if missing in restarts [mesh] @@ -44,170 +41,155 @@ Bxy = sqrt(Bpxy^2 + Btxy^2) bxcvz = 1./Rxy^2 # Curvature -paralleltransform = shifted - -symmetricGlobalX = true - -shiftangle_from_zshift=true +[mesh:paralleltransform] +type = shifted ################################################## # derivative methods -[ddx] +[mesh:ddx] first = C2 second = C2 upwind = W3 -[ddy] +[mesh:ddy] first = C2 second = C2 upwind = W3 -[ddz] +[mesh:ddz] -first = C2 -second = C2 +first = FFT +second = FFT upwind = W3 ################################################### # Time-integration solver [solver] - # Note: If evolving neutrals, need preconditioning type = cvode use_precon = true -ATOL = 1.0e-10 # absolute tolerance -RTOL = 1.0e-5 # relative tolerance +atol = 1e-10 # absolute tolerance +rtol = 1e-05 # relative tolerance mxstep = 1000000 # Maximum internal steps per output cvode_max_order = 2 cvode_stability_limit_detection = true -[laplace] # This is used for Delp2 operator -all_terms = true -nonuniform=true - ################################################## # Electrostatic potential solver [phiSolver] -inner_boundary_flags = 16 -outer_boundary_flags = 16 # INVERT_SET (2.8 * Te) - -all_terms = true -nonuniform=true # NOTE: Necessary to avoid numerical instability +inner_boundary_flags = 16 # INVERT_SET, setting inner boundary +outer_boundary_flags = 16 # INVERT_SET, setting outer boundary [laplacexy] # 2D solver in X-Y -pctype = sor # Preconditioner +pctype = hypre # Preconditioner atol = 1e-12 -rtol = 1e-8 +rtol = 1e-08 +# Note: These boundary flags should be consistent with phiSolver core_bndry_dirichlet = false pf_bndry_dirichlet = true y_bndry_dirichlet = false -include_y_derivs = true +include_y_derivs = true # Include the y components in the inversion [aparSolver] inner_boundary_flags = 0 outer_boundary_flags = 0 all_terms = true -nonuniform=true +nonuniform = true -# general settings for the model +################################################### +# Model settings [Hermes] ####################### # Output variables output_ddt = false # Output time derivatives -verbose = true # Output additional fields - -####################### -# radial buffers -radial_buffers = false -radial_inner_width = 4 -radial_outer_width = 4 -radial_buffer_D = 1.0 - -resistivity_boundary = 1e-2 -resistivity_boundary_width = 4 +verbose = true # Output more derived quantities ####################### # Numerical dissipation vepsi_dissipation = true # Parallel dissipation on Ve-Vi -vort_dissipation = true -hyperpar = 0.1 +vort_dissipation = false # Parallel dissipation of vorticity oscillations +phi_dissipation = true # Parallel dissipation of potential oscillations + +#x_hyper_viscos = 0.1 +z_hyper_viscos = 0.001 +ne_hyper_z = 0.001 +pe_hyper_z = 0.001 + +radial_buffers = true # Turn on dissipation close to boundaries +radial_inner_width = 4 # Number of cells on inner boundary +radial_outer_width = 4 # Number of cell on outer boundary. < 0 -> no buffer +radial_buffer_D = 1 # Damping coefficient in boundary + +####################### # Flux limiters kappa_limit_alpha = 0.2 # SOLPS style heat flux limiter -eta_limit_alpha = 0.5 # SOLPS style viscosity limiter +eta_limit_alpha = 0.5 # SOLPS style viscosity limiter ####################### # Electric field and Ohm's law electromagnetic = true # Electromagnetic? Otherwise electrostatic -FiniteElMass = true # Finite electron mass? +FiniteElMass = true # Finite electron mass? # Electrostatic potential -split_n0 = false # Solve n=0 separately? -split_n0_psi = false +split_n0 = false # Solve n=0 electrostatic potential separately? +split_n0_psi = false # Solve n=0 Apar component separately? -# NOTE: all currents switched off for fluid run -j_diamag = true # Diamagnetic current: Vort <-> Pe -j_par = true # Parallel current: Vort <-> Psi +phi_boundary_relax = true # Relax the boundary towards zero gradient +phi_boundary_timescale = 1e-05 # Timescale for relaxation [seconds] -pe_par = true # Parallel pressure gradient: Pe <-> Psi -resistivity = true # Resistivity: Psi -> Pe -thermal_flux = true -thermal_force = true -electron_viscosity = true -ion_viscosity = true # Ion parallel viscosity -thermal_conduction = true +####################### +# Terms included in the model -frecycle = 0.99 # Neutral gas recycling fraction +# NOTE: all currents switched off for fluid run +j_diamag = true # Diamagnetic current: Vort <-> Pe +j_par = true # Parallel current: Vort <-> Psi -carbon_fraction = 0.0 +pe_par = true # Parallel pressure gradient: Pe <-> Psi +resistivity = true # Resistivity: Psi -> Pe +thermal_flux = true +thermal_force = true +electron_viscosity = true +ion_viscosity = false # Viscous terms in vorticity equation +ion_viscosity_par = true # Parallel diffusion of ion momentum +thermal_conduction = true # Parallel heat conduction +poloidal_flows = false # Include ExB flows in the X-Y plane? +ion_velocity = true # Evolve ion parallel velocity? +boussinesq = true # Use Boussinesq approximation +ion_neutral = false # Include ion-neutral collisions in tau_i? +neutral_friction = false # Friction between plasma and neutrals in vorticity equation -excitation = true # Hydrogen neutral excitation radiation +frecycle = 0.99 # Neutral recycling fraction at targets -## Settings for 2D parallel closures -sinks = false -sink_invlpar = 0.2 # 5m parallel connection length -sheath_closure = false -drift_wave = false +carbon_fraction = -1 # Fixed fraction carbon impurity. < 0 means none ####################### # Transport coefficients -classical_diffusion = true # Collisional diffusion - -anomalous_D = -1 # Anomalous density diffusion [m^2/s] -anomalous_chi = -1 # Anomalous thermal diffusion [m^2/s] -anomalous_nu = -1 # Anomalous viscosity - -poloidal_flows = false - -magnetic_drift = true -ion_velocity = true - -ion_neutral = 0.0 -neutral_friction = false # Friction between plasma and neutrals +classical_diffusion = false # Collisional diffusion -boussinesq = true # Use Boussinesq approximation +anomalous_D = -1 # Anomalous density diffusion [m^2/s] +anomalous_chi = -1 # Anomalous thermal diffusion [m^2/s] +anomalous_nu = -1 # Anomalous viscosity [m^2/s] +####################### # Radial boundary fluxes -ne_bndry_flux = true -pe_bndry_flux = true -vort_bndry_flux = false - -ramp_mesh = false -ramp_timescale = 1e4 +ne_bndry_flux = true # Allow radial flows of particles through boundary +pe_bndry_flux = true # Allow radial flows of energy through boundary +vort_bndry_flux = false # Allow radial flows of vorticity through boundary ####################### # Plasma sheath @@ -216,38 +198,36 @@ sheath_yup = true sheath_ydown = false sheath_gamma_e = 4 # Electron sheath heat transmission sheath_gamma_i = 2.5 # Ion sheath heat transmission -neutral_gamma = 0.0 - -startprofiles = false +####################### +# Sources core_sources = false # Only sources in the core adapt_source = false # Feedback on profiles (PI controller) energy_source = true # Constant power per particle (at fixed x). False -> Constant power per volume -source_p = 1e-2 # Rate based on current error (1/timescale) -source_i = 1e-6 # Integral of error (1/timescale^2) +source_p = 0.01 # Rate based on current error (1/timescale) +source_i = 1e-06 # Integral of error (1/timescale^2) source_vary_g11 = false # Vary source in poloidal angle to better match radial transport -staggered = false - -# Normalisation parameters - +####################### # Normalisation factors -Nnorm = 1e20 +Nnorm = 1e+20 Tnorm = 100 Bnorm = 1 -AA = 1 # Atomic mass. 1 = Hydrogen, 2 = Deuterium +AA = 2 # Atomic mass. 1 = Hydrogen, 2 = Deuterium +####################### +# Neutrals [neutral] type = none # Neutral model: none, diffusion2d, recycling, fullvelocity, mixed viscosity = 1 # Dynamic viscosity bulk = 0 # Bulk (volume) viscosity -conduction = 1 -neutral_gamma = 0.0 +conduction = 1 +neutral_gamma = 0.0 # Sheath heat transmission -nn_floor = 1e-2 # Floor applied when calculating Vn = NVn / Nn -low_n_equilibriate = -1e-4 # If n < nn_floor, evolve Tn and Vn towards plasma values +nn_floor = 0.01 # Floor applied when calculating Vn = NVn / Nn +low_n_equilibriate = -0.0001 # If n < nn_floor, evolve Tn and Vn towards plasma values -[All] +[all] scale = 0.0 bndry_all = neumann_o2 @@ -265,7 +245,7 @@ function = 0 [VePsi] # Ve + 0.5*beta_e*mi_me*psi bndry_core = zerolaplace -bndry_pf = dirichlet_o2 +bndry_pf = dirichlet_o2 bndry_xout = dirichlet_o2 [Pe] # Electron pressure @@ -280,23 +260,19 @@ function = 0.01 source = 1e4 * H(x-0.5) * exp(-(x-0.5)/0.2) * H(0.2 - y/(2π)) -[Ve] - [phi] # Radial boundaries determined by Laplacian inversion -bndry_xin = none -bndry_xout = none - -bndry_all = dirichlet_o2 - - +bndry_all = none +bndry_yup = free_o3 +bndry_ydown = neumann -[Nn] +[Nn] # Neutral density -scale = 5e-2 +scale = 0.05 function = 1 +bndry_all = neumann_o2 -[Vn] +[Vn] # Neutral parallel velocity scale = 0 function = 0 @@ -318,13 +294,12 @@ scale = 0 function = 0 bndry_all = dirichlet_o2 -[Pn] +[Pn] # Neutral pressure -scale = 1e-5 +scale = 1e-05 function = 1 bndry_all = neumann_o2 [NVn] bndry_all = dirichlet_o2 - diff --git a/tokamak/1-no-currents/BOUT.inp b/tokamak/1-no-currents/BOUT.inp new file mode 100644 index 0000000..712ef41 --- /dev/null +++ b/tokamak/1-no-currents/BOUT.inp @@ -0,0 +1,203 @@ + +NOUT = 100 +TIMESTEP = 1000 + +MZ = 1 # number of points in z direction (2^n + 1) +ZPERIOD = 1 # Number of times domain repeats in toroidal angle + +grid = tokamak-example.nc + +[mesh] +paralleltransform = shifted + +extrapolate_y = false # Use metrics constant into y guard cells + +################################################## +# derivative methods + +[mesh:ddx] + +first = C2 +second = C2 +upwind = W3 + +[mesh:ddy] + +first = C2 +second = C2 +upwind = W3 + +[mesh:ddz] + +first = FFT +second = FFT +upwind = W3 + +################################################### +# Time-integration solver + +[solver] +# Note: If evolving neutrals, need preconditioning +type = cvode +use_precon = true + +maxl = 10 + +ATOL = 1.0e-10 # absolute tolerance +RTOL = 1.0e-5 # relative tolerance +mxstep = 1000000 # Maximum internal steps per output + + +################################################### +# Model settings + +[Hermes] + +####################### +# Output variables +output_ddt = false # Output time derivatives +verbose = true # Output more derived quantities + +####################### +# Numerical dissipation + +radial_buffers = true # Turn on dissipation close to boundaries +radial_inner_width = 4 # Number of cells on inner boundary +radial_outer_width = -1 # Number of cell on outer boundary. < 0 -> no buffer +radial_buffer_D = 1 # Damping coefficient in boundary + +####################### +# Flux limiters +kappa_limit_alpha = 0.2 # SOLPS style heat flux limiter +eta_limit_alpha = 0.5 # SOLPS style viscosity limiter + +####################### +# Terms included in the model + +# NOTE: all currents switched off for fluid run +j_diamag = false # Diamagnetic current +j_par = false # Parallel current + +pe_par = true # Parallel pressure gradient terms +ion_viscosity = true # Parallel diffusion of ion momentum +thermal_conduction = true # Parallel heat conduction +ion_velocity = true # Evolve ion parallel velocity? +ion_neutral = false # Include ion-neutral collisions in tau_i? + +frecycle = 0.99 # Neutral recycling fraction at targets + +carbon_fraction = 0.01 # Fixed fraction carbon impurity + +####################### +# Transport coefficients +classical_diffusion = true # Collisional diffusion + +anomalous_D = 1 # Anomalous density diffusion [m^2/s] +anomalous_chi = 0.5 # Anomalous thermal diffusion [m^2/s] +anomalous_nu = 0.01 # Anomalous viscosity [m^2/s] + +####################### +# Radial boundary fluxes +ne_bndry_flux = true # Allow radial flows of particles through boundary +pe_bndry_flux = true # Allow radial flows of energy through boundary + +####################### +# Plasma sheath +sheath_model = 2 # 0 = Bohn, 1 = Loizu, 2 = Bohm + free density +sheath_yup = true +sheath_ydown = true +sheath_gamma_e = 4 # Electron sheath heat transmission +sheath_gamma_i = 2.5 # Ion sheath heat transmission + +####################### +# Sources +core_sources = true # Only sources in the core +adapt_source_p = false # Feedback on pressure profiles from grid file? (PI controller) +adapt_source_n = true # Feedback on density profile from grid file? +energy_source = true # Constant power per particle (at fixed x). False -> Constant power per volume +source_p = 1e-2 # Rate based on current error (1/timescale) +source_i = 1e-6 # Integral of error (1/timescale^2) +source_vary_g11 = true # Vary source in poloidal angle to better match radial transport + +####################### +# Normalisation factors +Nnorm = 1e20 +Tnorm = 100 +Bnorm = 1 +AA = 2 # Atomic mass. 1 = Hydrogen, 2 = Deuterium + +####################### +# Neutrals +[neutral] +type = mixed +viscosity = 1 # Dynamic viscosity +bulk = 0 # Bulk (volume) viscosity +conduction = 1 +neutral_gamma = 0.0 # Sheath heat transmission + +nn_floor = 1e-2 # Floor applied when calculating Vn = NVn / Nn +low_n_equilibriate = -1e-4 # If n < nn_floor, evolve Tn and Vn towards plasma values + +[All] +scale = 0.0 + +bndry_all = neumann_o2 +bndry_xin = neumann_o2 +bndry_xout = neumann_o2 + +[Ne] # Electron density +scale = 1 +function = 0.0 # Starting value, added to grid profiles + +source = 0.0 # Starting particle source + +[Pe] # Electron pressure +scale = 1 +function = 0.0 # Starting value, added to grid profiles + +source = 1e4 * exp(-(32*x^2)) # Source of thermal energy into electrons + +[Pi] +scale = 1 +function = 0.0 # Starting value, added to grid profiles + +source = 1e4 * exp(-(32*x^2)) # Source of thermal energy into ions + +[Nn] # Neutral density + +scale = 5e-2 +function = 1 +bndry_xin = dirichlet_o2(1e-2) +bndry_xout = dirichlet_o2(1e-2) + +[Vn] # Neutral parallel velocity + +scale = 0 +function = 0 + +bndry_all = dirichlet_o2 + +[Vn_x] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_y] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_z] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Pn] # Neutral pressure + +scale = 1e-5 +function = 1 + +bndry_all = neumann_o2 + +[NVn] # Neutral parallel momentum +bndry_all = dirichlet_o2 diff --git a/tokamak/2-currents/BOUT.inp b/tokamak/2-currents/BOUT.inp new file mode 100644 index 0000000..3e5caa9 --- /dev/null +++ b/tokamak/2-currents/BOUT.inp @@ -0,0 +1,273 @@ + +NOUT = 100 +TIMESTEP = 200 + +MZ = 1 # number of points in z direction (2^n + 1) +ZPERIOD = 1 # Number of times domain repeats in toroidal angle + +grid = tokamak-example.nc + +[restart] +init_missing = true # Set missing variables to 0 if missing in restarts + +[mesh] +paralleltransform = shifted + +extrapolate_y = false # Use metrics constant into y guard cells + +################################################## +# derivative methods + +[mesh:ddx] + +first = C2 +second = C2 +upwind = W3 + +[mesh:ddy] + +first = C2 +second = C2 +upwind = W3 + +[mesh:ddz] + +first = FFT +second = FFT +upwind = W3 + +################################################### +# Time-integration solver + +[solver] +# Note: If evolving neutrals, need preconditioning +type = cvode +use_precon = true + +maxl = 10 + +ATOL = 1.0e-10 # absolute tolerance +RTOL = 1.0e-5 # relative tolerance +mxstep = 1000000 # Maximum internal steps per output + +################################################## +# Electrostatic potential solver + +[phiSolver] +inner_boundary_flags = 0 # Zero-gradient inner boundary +outer_boundary_flags = 16 # INVERT_SET, setting outer boundary + +[laplacexy] # 2D solver in X-Y +pctype = hypre # Preconditioner + +atol = 1e-12 +rtol = 1e-8 + +# Note: These boundary flags should be consistent with phiSolver +core_bndry_dirichlet = false +pf_bndry_dirichlet = true +y_bndry_dirichlet = false + +include_y_derivs = true # Include the y components in the inversion + +[aparSolver] +inner_boundary_flags = 0 +outer_boundary_flags = 0 + +all_terms = true +nonuniform=true + +################################################### +# Model settings + +[Hermes] + +####################### +# Output variables +output_ddt = false # Output time derivatives +verbose = true # Output more derived quantities + +####################### +# Numerical dissipation + +vepsi_dissipation = true # Parallel dissipation on Ve-Vi +vort_dissipation = false # Parallel dissipation of vorticity oscillations +phi_dissipation = true # Parallel dissipation of potential oscillations + +radial_buffers = true # Turn on dissipation close to boundaries +radial_inner_width = 4 # Number of cells on inner boundary +radial_outer_width = -1 # Number of cell on outer boundary. < 0 -> no buffer +radial_buffer_D = 1 # Damping coefficient in boundary + +####################### +# Flux limiters +kappa_limit_alpha = 0.2 # SOLPS style heat flux limiter +eta_limit_alpha = 0.5 # SOLPS style viscosity limiter + +####################### +# Electric field and Ohm's law +electromagnetic = true # Electromagnetic? Otherwise electrostatic +FiniteElMass = true # Finite electron mass? + +# Electrostatic potential +split_n0 = false # Solve n=0 electrostatic potential separately? +split_n0_psi = false # Solve n=0 Apar component separately? + +phi_boundary_relax = true # Relax the boundary towards zero gradient +phi_boundary_timescale = 1e-5 # Timescale for relaxation [seconds] + +####################### +# Terms included in the model + +# NOTE: all currents switched off for fluid run +j_diamag = true +j_par = true + +pe_par = true # Parallel pressure gradient: Pe <-> Psi +resistivity = true # Resistivity: Psi -> Pe +thermal_flux = true +thermal_force = true +electron_viscosity = true +ion_viscosity = false # Viscous terms in vorticity equation +ion_viscosity_par = true # Parallel diffusion of ion momentum +thermal_conduction = true # Parallel heat conduction +poloidal_flows = false # Include ExB flows in the X-Y plane? +ion_velocity = true # Evolve ion parallel velocity? +boussinesq = true # Use Boussinesq approximation +ion_neutral = false # Include ion-neutral collisions in tau_i? +neutral_friction = false # Friction between plasma and neutrals in vorticity equation + +frecycle = 0.99 # Neutral recycling fraction at targets + +carbon_fraction = 0.01 # Fixed fraction carbon impurity + +####################### +# Transport coefficients +classical_diffusion = true # Collisional diffusion + +anomalous_D = 1 # Anomalous density diffusion [m^2/s] +anomalous_chi = 0.5 # Anomalous thermal diffusion [m^2/s] +anomalous_nu = 0.01 # Anomalous viscosity [m^2/s] + +####################### +# Radial boundary fluxes +ne_bndry_flux = true # Allow radial flows of particles through boundary +pe_bndry_flux = true # Allow radial flows of energy through boundary +vort_bndry_flux = false # Allow radial flows of vorticity through boundary + +####################### +# Plasma sheath +sheath_model = 2 # 0 = Bohn, 1 = Loizu, 2 = Bohm + free density +sheath_yup = true +sheath_ydown = true +sheath_gamma_e = 4 # Electron sheath heat transmission +sheath_gamma_i = 2.5 # Ion sheath heat transmission + +####################### +# Sources +core_sources = true # Only sources in the core +adapt_source_p = false # Feedback on pressure profiles from grid file? (PI controller) +adapt_source_n = true # Feedback on density profile from grid file? +energy_source = true # Constant power per particle (at fixed x). False -> Constant power per volume +source_p = 1e-2 # Rate based on current error (1/timescale) +source_i = 1e-6 # Integral of error (1/timescale^2) +source_vary_g11 = true # Vary source in poloidal angle to better match radial transport + +####################### +# Normalisation factors +Nnorm = 1e20 +Tnorm = 100 +Bnorm = 1 +AA = 2 # Atomic mass. 1 = Hydrogen, 2 = Deuterium + +####################### +# Neutrals +[neutral] +type = mixed +viscosity = 1 # Dynamic viscosity +bulk = 0 # Bulk (volume) viscosity +conduction = 1 +neutral_gamma = 0.0 # Sheath heat transmission + +nn_floor = 1e-2 # Floor applied when calculating Vn = NVn / Nn +low_n_equilibriate = -1e-4 # If n < nn_floor, evolve Tn and Vn towards plasma values + +[All] +scale = 0.0 + +bndry_all = neumann_o2 +bndry_xin = neumann_o2 +bndry_xout = neumann_o2 + +[Ne] # Electron density +scale = 1 +function = 0.0 # Starting value, added to grid profiles + +source = 0.0 # Starting particle source + +[Vort] +function = 0 + +[VePsi] # Ve + 0.5*beta_e*mi_me*psi +bndry_core = zerolaplace +bndry_pf = dirichlet_o2 +bndry_xout = dirichlet_o2 + +[Pe] # Electron pressure +scale = 1 +function = 0.0 # Starting value, added to grid profiles + +source = 1e4 * exp(-(32*x^2)) # Source of thermal energy into electrons + +[Pi] +scale = 1 +function = 0.0 # Starting value, added to grid profiles + +source = 1e4 * exp(-(32*x^2)) # Source of thermal energy into ions + +[Ve] + +[phi] +# Radial boundaries determined by Laplacian inversion +bndry_all = none +bndry_yup = free_o3 +bndry_ydown = free_o3 + +[Nn] # Neutral density + +scale = 5e-2 +function = 1 +bndry_xin = dirichlet_o2(1e-2) +bndry_xout = dirichlet_o2(1e-2) + +[Vn] # Neutral parallel velocity + +scale = 0 +function = 0 + +bndry_all = dirichlet_o2 + +[Vn_x] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_y] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_z] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Pn] # Neutral pressure + +scale = 1e-5 +function = 1 + +bndry_all = neumann_o2 + +[NVn] # Neutral parallel momentum +bndry_all = dirichlet_o2 diff --git a/tokamak/README.md b/tokamak/README.md new file mode 100644 index 0000000..514703a --- /dev/null +++ b/tokamak/README.md @@ -0,0 +1,89 @@ +Tokamak example +=============== + +Input grid +---------- + +To run these examples, you will need a grid file for a tokamak. +These can be generated by [Hypnotoad](https://github.com/boutproject/hypnotoad). + +The BOUT.inp input files assume that this grid file is in the same +directory as this README, and is called `tokamak-grid.nc` + +Plasma profiles +--------------- + +To run a simulation, the density and temperature profiles should be +set in the grid file. These profiles are used for: +- Initial conditions. By default the `Ne0`, `Te0` and `Ti0` profiles are used + to add to the starting `Ne`, `Pe` and `Pi` fields. This can be disabled by + setting `startprofiles` to `false` in the `BOUT.inp` file. +- Feedback control of the particle source, if `adapt_source_n` is `true`. This uses + the difference between the electron density and the value in the grid file. +- Feedback control of the electron and ion power sources, if `adapt_source_n` is `true`. + These use the difference between the electron/ion pressure and the grid file profiles + to set the electron/ion input power. + +To modify a grid file, specifying the separatrix temperature and density: + + $ ./setprofiles thegridfile.nc Tesep Nesep + +where `Tesep` is in eV and `Nesep` is in m^-3. e.g + + $ ./setprofiles thegridfile.nc 100 1e19 + +This will modify in place `thegridfile.nc`, keeping the pressure gradient +unchanged. + +Curvature +--------- + +The magnetic curvature is used in calculating the drift +when diamagnetic current is enabled. To ensure that the flux-surface +average radial flux is zero when pressure is constant on a flux +surface, the quantity J * Curl(b/B) must be zero. If this is not the case, +then large spurious radial currents and flows can result. + +A python script is included, `adjust_curvature.py` which will modify +the radial component of the curvature such that J * Curl(b/B) sums to +zero over closed flux surfaces. Note that this doesn't currently ensure that +Div( Curl(b/B) ) is zero, as the poloidal component of Curl(b/B) is not +modified. To use the script + + $ ./adjust_curvature.py thegridfile.nc + +This will modify `thegridfile.nc`. + +Simulations +----------- + +First run the transport case. This includes neutrals, cross-field +diffusion, but no currents or drifts. + + $ mpirun -np 16 ./hermes-2 -d 1-no-currents + +This will run for approx. 30 minutes, for a simulation time of 1e5/wci, +sufficient to reach steady state. + +Then turn on currents and diamagnetic drift. Copy the `BOUT.restart.*` files +from `1-no-currents` into `2-currents` and restart: + + $ cp 1-no-currents/BOUT.restart.* 2-currents/ + $ mpirun -np 16 ./hermes-2 -d 2-currents restart + +Analysis +-------- + +There are some python scripts to make plots and perform analysis of the volume +integrated power sources, sinks and transfer terms. Most of these take two +arguments: The grid file, and the directory containing the data e.g. + + $ ./make2dplots.py tokamak-example.nc 1-no-currents + + +- `make2dplots.py` Make 2D (poloidal cross-section) plots of the density and + temperature; radiation; and energy transfer +- `plot_target_profiles.py` Plot density, temperature and heat flux along + the inner and outer targets +- `volume_integrals.py` Integrate the power sources, sinks, and transfer terms + over the core and SOL regions of the grid. diff --git a/tokamak/adjust_curvature.py b/tokamak/adjust_curvature.py new file mode 100755 index 0000000..06b540f --- /dev/null +++ b/tokamak/adjust_curvature.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 +# +# Check the poloidal integral of the radial magnetic drift, Curl(b/B) +# +# The sum over the core poloidal points of J * Curlb_Bx should be zero + +import sys + +if len(sys.argv) != 2: + print(f"Usage: {sys.argv[0]} gridfile.nc") + print("Note: The input grid file is modified") + sys.exit(1) + +filename = sys.argv[1] + +from boututils.datafile import DataFile +import numpy as np + +# Read values from current grid file +with DataFile(filename) as g: + J = g["J"] + B = g["Bxy"] + + # Curvature components + bxcvx = g["bxcvx"] + bxcvz = g["bxcvz"] + + sinty = g["sinty"] # Integrated magnetic shear + + # Curl(b / B) + curlb_Bx = 2 * bxcvx / B + + # Poloidal branch-cut indices + j11 = g["jyseps1_1"] + j12 = g["jyseps1_2"] + j21 = g["jyseps2_1"] + j22 = g["jyseps2_2"] + + # Get the separatrix index + ix1 = g["ixseps1"] + ix2 = g["ixseps2"] + ixcore = np.min([ix1, ix2]) + + # The coordinate system depends on the transform + try: + parallel_transform = g["parallel_transform"] + except KeyError: + parallel_transform = "identity" + +# Calculate J * Curl(b/B) in the core region +Jc = J * curlb_Bx +Jc_core = np.concatenate( + (Jc[:, (j11 + 1) : (j21 + 1)], Jc[:, (j12 + 1) : (j22 + 1)]), axis=1 +) + +# Sum over y (poloidal) index in the core +Jc_sum = np.sum(Jc_core, axis=1) + +# Rescale curvature +# Keep places where curvature = 0, so multiply positive and negative values +# by factors. +# Positive values are multiplied by "factor"; negative values are divided +# +# Curvature = positive_curvature * factor + negative_curvature / factor +# + + +# Sum over just the positive part of the curvature +Jc_positive_sum = np.sum(np.clip(Jc_core, 0.0, None), axis=1) + +factor = np.sqrt(1 - Jc_sum / Jc_positive_sum) + +# Keep the factor constant outside the separatrix +# so that curvature is still continuous +factor_edge = factor[ixcore - 1] # Last closed flux surface + +bxcvx2 = ( + np.clip(bxcvx, 0.0, None) * factor_edge + + np.clip(bxcvx, None, 0.0) / factor_edge # Positive parts +) # Negative parts + +# Go through core flux surfaces +for i in range(ixcore): + # Inner core + bxcvx2[i, (j11 + 1) : (j21 + 1)] = ( + np.clip(bxcvx[i, (j11 + 1) : (j21 + 1)], 0.0, None) * factor[i] + + np.clip(bxcvx[i, (j11 + 1) : (j21 + 1)], None, 0.0) / factor[i] + ) + + # Outer core + bxcvx2[i, (j12 + 1) : (j22 + 1)] = ( + np.clip(bxcvx[i, (j12 + 1) : (j22 + 1)], 0.0, None) * factor[i] + + np.clip(bxcvx[i, (j12 + 1) : (j22 + 1)], None, 0.0) / factor[i] + ) + + +# Check that the sum of the curvature sums to zero over poloidal angle + +Jc2 = J * 2 * bxcvx2 / B +Jc2_core = np.concatenate( + (Jc2[:, (j11 + 1) : (j21 + 1)], Jc2[:, (j12 + 1) : (j22 + 1)]), axis=1 +) +Jc2_sum = np.sum(Jc2_core, axis=1) + +import matplotlib.pyplot as plt + +plt.plot(Jc_sum, label="Input ") +plt.plot(Jc2_sum, label="Adjusted") +plt.xlabel("Radial index") +plt.ylabel(r"$\sum_{y,core} J \nabla\times(\mathbf{b}/B)$") +plt.legend() +plt.show() + +with DataFile(filename, write=True) as outf: + # Update curvature + outf["bxcvx"] = bxcvx2 + + if parallel_transform == "identity": + # Need to correct the z component + outf["bxcvz"] = bxcvz + sinty * (bxcvx - bxcvx2) diff --git a/tokamak/hermes-2 b/tokamak/hermes-2 new file mode 120000 index 0000000..d88e98a --- /dev/null +++ b/tokamak/hermes-2 @@ -0,0 +1 @@ +../hermes-2 \ No newline at end of file diff --git a/tokamak/make2dplots.py b/tokamak/make2dplots.py new file mode 100755 index 0000000..493faa3 --- /dev/null +++ b/tokamak/make2dplots.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python +# +# Makes 2D plots: +# - ne_nn_te_2d.png Electron density, neutral density and electron temperature +# - radiation2d.png Impurity and hydrogenic radiation +# - transfer2d.png Ion <-> Electron and Ion <-> Neutral energy exchange +# + +import argparse + +parser = argparse.ArgumentParser(description="Make 2D plots of Hermes-2 data") +parser.add_argument("gridfile", type=str, help="The grid file used") +parser.add_argument("datapath", type=str, help="The path to the data (BOUT.dmp files)") + +args = parser.parse_args() +gridfile = args.gridfile +path = args.datapath + +from boutdata import collect +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib.ticker import LogFormatter + +from boututils.datafile import DataFile +from boutdata.griddata import gridcontourf +import numpy as np + +t_arr = collect("t_array", path=path) +tind = len(t_arr) - 2 + +Nnorm = collect("Nnorm", path=path) +Tnorm = collect("Tnorm", path=path) +wci = collect("Omega_ci", path=path) + +grid = DataFile(gridfile) + +ni = collect("Ne", path=path, tind=tind).squeeze() * Nnorm # m^-3 +nn = collect("Nn", path=path, tind=tind).squeeze() * Nnorm # m^-3 +te = collect("Telim", path=path, tind=tind).squeeze() * Tnorm # eV + +####################### + +fig, ax = plt.subplots(1, 3, figsize=(10, 8)) + +c = gridcontourf(grid, ni, ax=ax[0], show=False) +fig.colorbar(c, ax=ax[0]) +ax[0].set_title(r"Ion dens [m$^{3}$]") +ax[0].axes.get_xaxis().set_ticks([1, 2, 3, 4]) + +c = gridcontourf(grid, nn, ax=ax[1], show=False, ylabel=None) +fig.colorbar(c, ax=ax[1]) +ax[1].set_title(r"Atom dens [m$^{3}$]") +ax[1].axes.get_yaxis().set_ticks([]) +ax[1].axes.get_xaxis().set_ticks([1, 2, 3, 4]) + +c = gridcontourf(grid, te, ax=ax[2], log=True, nlevel=20, show=False, ylabel=None) +cbar = fig.colorbar(c, ax=ax[2]) +cbar.set_ticks([0.1, 0.5, 1.0, 5.0, 10.0, 50, 100, 200, 500]) +ax[2].set_title("Elec. temp [eV]") +ax[2].axes.get_yaxis().set_ticks([]) +ax[2].axes.get_xaxis().set_ticks([1, 2, 3, 4]) + +fig.tight_layout(pad=3.0) + +plt.savefig("ne_nn_te_2d.png") +plt.show() +plt.close("all") + +####################### + +rp = collect("Rp", path=path, tind=tind).squeeze() * Nnorm * Tnorm * 1.602e-19 * wci + +try: + rzrad = ( + collect("Rzrad", path=path, tind=tind).squeeze() + * Nnorm + * Tnorm + * 1.602e-19 + * wci + ) # W/m^3 +except: + rzrad = np.zeros(rp.shape) + +# Same scale for all plots +maxd = max(np.max(rzrad), np.max(rp)) +mind = min(np.min(rzrad), np.min(rp)) + +fig, ax = plt.subplots(2, 1, figsize=(10, 7), gridspec_kw={"height_ratios": [10, 1]}) +from mpl_toolkits.axes_grid1 import make_axes_locatable + +ne_ax = ax[0] +cbar_ax = ax[1] +divider = make_axes_locatable(ne_ax) +h_ax = divider.append_axes("right", size="100%", pad=0.01) + +c = gridcontourf(grid, rzrad, ax=ne_ax, show=False, mind=mind, maxd=maxd) +ne_ax.set_title(r"Impurity radiation [Wm$^{-3}$]") +ne_ax.axes.get_xaxis().set_ticks([1, 2, 3, 4]) +ne_ax.set_ylim(bottom=4.0) + +ne_ax.xaxis.grid(True) +ne_ax.yaxis.grid(True) + +c = gridcontourf(grid, rp, ax=h_ax, show=False, ylabel=None, mind=mind, maxd=maxd) +h_ax.set_title(r"H radiation [Wm$^{-3}$]") +h_ax.axes.get_yaxis().set_ticklabels([]) +h_ax.axes.get_xaxis().set_ticks([1, 2, 3, 4]) +h_ax.set_ylim(bottom=4.0) + +h_ax.xaxis.grid(True) +h_ax.yaxis.grid(True) + +cbar = fig.colorbar(c, cax=cbar_ax, orientation="horizontal") + +fig.tight_layout() + +plt.savefig("radiation2d.png") + +plt.show() +plt.close("all") + +############################### +# Power transfer + +qi = collect("Qi", path=path, tind=tind).squeeze() * Nnorm * Tnorm * 1.602e-19 * wci +wi = ( + collect("Wi", path=path, tind=tind).squeeze() * Nnorm * Tnorm * 1.602e-19 * wci +) # elec -> ion + +cmap = plt.cm.get_cmap("bwr") + +maxd = max([np.max(np.abs(qi)), np.max(np.abs(wi))]) * 1e-6 +mind = -maxd + +fig, ax = plt.subplots(2, 1, figsize=(10, 7), gridspec_kw={"height_ratios": [10, 1]}) +from mpl_toolkits.axes_grid1 import make_axes_locatable + +qi_ax = ax[0] +cbar_ax = ax[1] +divider = make_axes_locatable(qi_ax) +wi_ax = divider.append_axes("right", size="100%", pad=0.01) + +c = gridcontourf(grid, qi * 1e-6, ax=qi_ax, show=False, mind=mind, maxd=maxd, cmap=cmap) +qi_ax.set_title(r"Ion$\rightarrow$Neutral transfer [MWm$^{-3}$]") +qi_ax.axes.get_xaxis().set_ticks([1, 2, 3, 4]) +qi_ax.set_ylim(bottom=4.0) + +qi_ax.xaxis.grid(True) +qi_ax.yaxis.grid(True) + +c = gridcontourf( + grid, -wi * 1e-6, ax=wi_ax, show=False, ylabel=None, mind=mind, maxd=maxd, cmap=cmap +) +wi_ax.set_title(r"Ion$\rightarrow$Electron transfer [MWm$^{-3}$]") +wi_ax.axes.get_yaxis().set_ticklabels([]) +wi_ax.axes.get_xaxis().set_ticks([1, 2, 3, 4]) +wi_ax.set_ylim(bottom=4.0) + +wi_ax.xaxis.grid(True) +wi_ax.yaxis.grid(True) + +cbar = fig.colorbar(c, cax=cbar_ax, orientation="horizontal") + +fig.tight_layout() + +plt.savefig("transfer2d.png") + +plt.show() +plt.close("all") diff --git a/tokamak/plot_target_profiles.py b/tokamak/plot_target_profiles.py new file mode 100755 index 0000000..7d5b8ae --- /dev/null +++ b/tokamak/plot_target_profiles.py @@ -0,0 +1,209 @@ +#!/usr/bin/env python +# +# Make plots of the density, temperature and heat flux on the inner and outer targets +# Saves to outer_target_profiles and outer_target_profiles PDF and PNG files + +import argparse + +parser = argparse.ArgumentParser(description="Make plots of target profiles") +parser.add_argument("gridfile", type=str, help="The grid file used") +parser.add_argument("datapath", type=str, help="The path to the data (BOUT.dmp files)") + +args = parser.parse_args() +gridfile = args.gridfile +path = args.datapath + +# sheath heat transmission +gamma_e = 4.0 +gamma_i = 2.5 + +from boutdata import collect + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +from boututils.datafile import DataFile +import numpy as np + + +with DataFile(gridfile) as grid: + Rxy = grid["Rxy"] + Zxy = grid["Zxy"] + Bpxy = grid["Bpxy"] + Bxy = grid["Bxy"] + ixsep = grid["ixseps1"] + +J = collect("J", path=path) +g_22 = collect("g_22", path=path) +dx = collect("dx", path=path) +dy = collect("dy", path=path) +dz = 2 * np.pi + +dV = abs(J * dx * dy * dz) # Volume element +dS = dV / (dy * np.sqrt(g_22)) # Area element + +# Calculate distance along target + +# Outer target +dR = Rxy[3:-1, -1] - Rxy[2:-2, -1] +dZ = Zxy[3:-1, -1] - Zxy[2:-2, -1] + +dl = np.sqrt(dR ** 2 + dZ ** 2) +distance = np.cumsum(dl) + +# Inner target +dR = Rxy[3:-1, 0] - Rxy[2:-2, 0] +dZ = Zxy[3:-1, 0] - Zxy[2:-2, 0] + +dl_inner = np.sqrt(dR ** 2 + dZ ** 2) +distance_inner = np.cumsum(dl_inner) + +# Set to zero at separatrix, which is half-way between +# grid centers +distance -= 0.5 * (distance[ixsep] + distance[ixsep - 1]) +distance_inner -= 0.5 * (distance_inner[ixsep] + distance_inner[ixsep - 1]) + +t_arr = collect("t_array", path=path) +tind = len(t_arr) - 2 + +Nnorm = collect("Nnorm", path=path) +Tnorm = collect("Tnorm", path=path) +wci = collect("Omega_ci", path=path) +cs0 = collect("Cs0", path=path) +rho_s = collect("rho_s0", path=path) + +ni = collect("Ne", path=path, tind=tind, yguards=True).squeeze() * Nnorm # m^-3 +nn = collect("Nn", path=path, tind=tind, yguards=True).squeeze() * Nnorm # m^-3 +te = collect("Telim", path=path, tind=tind, yguards=True).squeeze() * Tnorm # eV +ti = collect("Ti", path=path, tind=tind, yguards=True).squeeze() * Tnorm # eV +pn = collect("Pn", path=path, tind=tind, yguards=True).squeeze() +tn = pn / (nn / Nnorm) * Tnorm # eV + +# Outer target +ni_targ = 0.5 * (ni[:, -3] + ni[:, -2])[2:-2] +nn_targ = 0.5 * (nn[:, -3] + nn[:, -2])[2:-2] +te_targ = 0.5 * (te[:, -3] + te[:, -2])[2:-2] +ti_targ = 0.5 * (ti[:, -3] + ti[:, -2])[2:-2] +tn_targ = 0.5 * (tn[:, -3] + tn[:, -2])[2:-2] + +cs = np.sqrt((ti_targ + te_targ) / Tnorm) * cs0 # m/s + +q_heat_e = gamma_e * 1.602e-19 * te_targ * ni_targ * cs # W/m^2 +q_heat_i = gamma_i * 1.602e-19 * ti_targ * ni_targ * cs # W/m^2 +q_pot = ni_targ * cs * 1.602e-19 * (13.6 + 2.2) # Potential energy + +# Inner target + +ni_targ_inner = 0.5 * (ni[:, 1] + ni[:, 2])[2:-2] +nn_targ_inner = 0.5 * (nn[:, 1] + nn[:, 2])[2:-2] +te_targ_inner = 0.5 * (te[:, 1] + te[:, 2])[2:-2] +ti_targ_inner = 0.5 * (ti[:, 1] + ti[:, 2])[2:-2] +tn_targ_inner = 0.5 * (tn[:, 1] + tn[:, 2])[2:-2] + +cs_inner = np.sqrt((ti_targ_inner + te_targ_inner) / Tnorm) * cs0 # m/s + +q_heat_e_inner = gamma_e * 1.602e-19 * te_targ_inner * ni_targ_inner * cs_inner # W/m^2 +q_heat_i_inner = gamma_i * 1.602e-19 * ti_targ_inner * ni_targ_inner * cs_inner # W/m^2 +q_pot_inner = ni_targ_inner * cs_inner * 1.602e-19 * (13.6 + 2.2) # Potential energy + +########### + +# Convert parallel power flux to perpendicular +qperp = (q_heat_e + q_heat_i + q_pot) * np.abs(Bpxy[2:-2, -1] / Bxy[2:-2, -1]) + +# Integrate to get total power +power_total = np.trapz(qperp * dl * 2 * np.pi * Rxy[2:-2, -1]) +power_total2 = np.sum( + (q_heat_e + q_heat_i + q_pot) * dS[2:-2, -1] * rho_s ** 2 +) # Cross-check + +print("Total power outer target power: {} MW".format(power_total * 1e-6)) +print("Total power outer target power (alt.): {} MW".format(power_total2 * 1e-6)) + + +qperp = (q_heat_e_inner + q_heat_i_inner + q_pot_inner) * np.abs( + Bpxy[2:-2, 0] / Bxy[2:-2, 0] +) +power_total = np.trapz(qperp * dl_inner * 2 * np.pi * Rxy[2:-2, 0]) +print("Total power inner target power: {} MW".format(power_total * 1e-6)) + +############ +# Outer target + +fig, ax = plt.subplots(1, 3, figsize=(10, 4)) + +ax[0].plot(distance, ni_targ, "-k", label="Ion") +ax[0].plot(distance, nn_targ, "--b", label="Atom") +ax[0].set_yscale("log") +ax[0].set_xlabel("Perpendicular distance [m]") +ax[0].axvline(0.0, linestyle="--", color="k") +ax[0].legend() +ax[0].set_ylabel(r"Density [m$^{-3}$]") + +ax[1].plot(distance, ti_targ, "-k", label="Ion") +ax[1].plot(distance, te_targ, "-r", label="Electron") +ax[1].plot(distance, tn_targ, "--b", label="Atom") +ax[1].set_ylabel(r"Temperature [eV]") +ax[1].set_xlabel("Perpendicular distance [m]") +ax[1].axvline(0.0, linestyle="--", color="k") +ax[1].set_ylim(bottom=0.0) +ax[1].legend() + +ax[2].plot(distance, q_heat_i * 1e-6, "-k", label="Ion") +ax[2].plot(distance, q_heat_e * 1e-6, "-r", label="Electron") +ax[2].plot(distance, q_pot * 1e-6, "-g", label="Potential") +ax[2].set_ylabel(r"Parallel heat flux [MW/m$^2$]") +ax[2].set_xlabel("Perpendicular distance [m]") +ax[2].set_ylim(bottom=0.0) +ax[2].axvline(0.0, linestyle="--", color="k") +ax[2].legend() + +fig.tight_layout() + +plt.savefig("outer_target_profiles.pdf") +plt.savefig("outer_target_profiles.png") +plt.close("all") + +############ +# Inner target + +fig, ax = plt.subplots(1, 3, figsize=(10, 4)) + +ax[0].plot(distance_inner, ni_targ_inner, "-k", label="Ion") +ax[0].plot(distance_inner, nn_targ_inner, "--b", label="Atom") +ax[0].set_yscale("log") +ax[0].set_xlabel("Perpendicular distance [m]") +ax[0].axvline(0.0, linestyle="--", color="k") +ax[0].legend() +ax[0].set_ylabel(r"Density [m$^{-3}$]") + +ax[1].plot(distance_inner, ti_targ_inner, "-k", label="Ion") +ax[1].plot(distance_inner, te_targ_inner, "-r", label="Electron") +ax[1].plot(distance_inner, tn_targ_inner, "--b", label="Atom") +ax[1].set_ylabel(r"Temperature [eV]") +ax[1].set_xlabel("Perpendicular distance [m]") +ax[1].axvline(0.0, linestyle="--", color="k") +ax[1].set_ylim(bottom=0.0) +ax[1].legend() + +ax[2].plot(distance_inner, q_heat_i_inner * 1e-6, "-k", label="Ion") +ax[2].plot(distance_inner, q_heat_e_inner * 1e-6, "-r", label="Electron") +ax[2].plot(distance_inner, q_pot_inner * 1e-6, "-g", label="Potential") +ax[2].set_ylabel(r"Parallel heat flux [MW/m$^2$]") +ax[2].set_xlabel("Perpendicular distance [m]") +ax[2].set_ylim(bottom=0.0) +ax[2].axvline(0.0, linestyle="--", color="k") +ax[2].legend() + +xticks = [-0.05, 0.0, 0.05] +for i in range(3): + ax[i].set_xticks(xticks) + + +fig.tight_layout() + +plt.savefig("inner_target_profiles.pdf") +plt.savefig("inner_target_profiles.png") +plt.close("all") diff --git a/tokamak/setprofiles.py b/tokamak/setprofiles.py new file mode 100755 index 0000000..2749812 --- /dev/null +++ b/tokamak/setprofiles.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 +# +# Change profiles in input grid +# - Keep pressure gradient the same, so Grad-Shafranov is unchanged +# - Set separatrix density and temperature + +import argparse + +parser = argparse.ArgumentParser(description="Set profiles in a BOUT++ grid") +parser.add_argument("gridfile", type=str, help="The file to modify") +parser.add_argument("Tsep", type=float, help="Separatrix temperature [eV]") +parser.add_argument("Nsep", type=float, help="Separatrix density [m^-3]") + +args = parser.parse_args() +filename = args.gridfile +Tsep = args.Tsep +Nesep = args.Nsep + +from boututils.datafile import DataFile +import numpy as np + +qe = 1.602e-19 + +with DataFile(filename) as f: + P = 0.0 + if "pressure" in f.keys(): + P = f["pressure"] # In Pascals + # Check that P's not zero. This is because some grid files + # may have pressure = 0. + + if np.amax(P) < 1e-3: + print("WARNING: input grid has no pressure, or pressure is zero") + Ne = f["Ni0"] * 1e20 # Note: Grid file in 1e20 m^-3 + Te = f["Te0"] + Ti = f["Ti0"] + P = qe * Ne * (Te + Ti) + + ixseps = f["ixseps1"] + # Y index around the midplane. Doesn't matter as long as it's in the core region + yind = int(0.5 * (f["jyseps1_2"] + f["jyseps2_2"] + 1)) + +# Subtract pressure, to get required value at separatrix +# Clip to impose minimum of zero +P = np.clip(P - (P[ixseps, yind] - qe * Nesep * 2.0 * Tsep), 1e-3, None) + +shape = np.sqrt(P) / np.sqrt(P[ixseps, yind]) # 1 at the separatrix + +Te = shape * Tsep +Ti = Te +Ne = shape * (Nesep * 1e-20) + +with DataFile(filename, write=True, create=False) as f: + f["Ni0"] = Ne + f["Ne0"] = Ne + f["Te0"] = Te + f["Ti0"] = Ti diff --git a/tokamak/volume_integrals.py b/tokamak/volume_integrals.py new file mode 100755 index 0000000..453413e --- /dev/null +++ b/tokamak/volume_integrals.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python +# +# Calculate and print volume integrated quantities from Hermes-2 simulation + +import argparse + +parser = argparse.ArgumentParser( + description="Calculate volume integrated quantities from Hermes-2 simulation" +) +parser.add_argument("gridfile", type=str, help="The grid file used") +parser.add_argument("datapath", type=str, help="The path to the data (BOUT.dmp files)") + +args = parser.parse_args() +gridfile = args.gridfile +path = args.datapath + +from boutdata import collect +from boututils.datafile import DataFile +import numpy as np + +t_arr = collect("t_array", path=path) +tind = len(t_arr) - 2 + +Nnorm = collect("Nnorm", path=path) +Tnorm = collect("Tnorm", path=path) +wci = collect("Omega_ci", path=path) +rho_s = collect("rho_s0", path=path) + +power_norm = Nnorm * Tnorm * 1.602e-19 * rho_s ** 3 * wci # Power [W] +pdens_norm = Nnorm * Tnorm * 1.602e-19 * wci # Power density [W/m^3] + +J = collect("J", path=path)[2:-2, :] # Note: Remove X guard cells +dx = collect("dx", path=path)[2:-2, :] +dy = collect("dy", path=path)[2:-2, :] +dz = 2 * np.pi + +# Volume element +dV = np.abs(J * dx * dy * dz) + +# Create a mask which is 1 in the core +# Excludes X guard cells +def core_mask(grid): + ixseps = min([grid["ixseps1"], grid["ixseps2"]]) # Innermost separatrix + j11 = grid["jyseps1_1"] + j12 = grid["jyseps1_2"] + j21 = grid["jyseps2_1"] + j22 = grid["jyseps2_2"] + + mask = np.zeros((grid["nx"] - 4, grid["ny"])) + # Inner core + mask[: ixseps - 2, (j11 + 1) : (j21 + 1)] = 1.0 + # Outer core + mask[: ixseps - 2, (j12 + 1) : (j22 + 1)] = 1.0 + return mask + + +with DataFile(gridfile) as grid: + ixseps = grid["ixseps1"] + mask = core_mask(grid) + +Spe = collect("Spe", path=path, tind=tind).squeeze()[2:-2, :] +Spi = collect("Spi", path=path, tind=tind).squeeze()[2:-2, :] +Sn = collect("Sn", path=path, tind=tind).squeeze()[2:-2, :] + +try: + Rzrad = collect("Rzrad", path=path, tind=tind).squeeze()[ + 2:-2, : + ] # Impurity radiation +except: + Rzrad = np.zeros(Sn.shape) +Rp = collect("Rp", path=path, tind=tind).squeeze()[ + 2:-2, : +] # Hydrogenic radiation + ionization cost +Qi = collect("Qi", path=path, tind=tind).squeeze()[2:-2, :] # Transfer ions -> neutrals +Wi = collect("Wi", path=path, tind=tind).squeeze()[ + 2:-2, : +] # Transfer electrons -> ions + +Spe_tot = np.sum(Spe * dV) +Spi_tot = np.sum(Spi * dV) +Sn_tot = np.sum(Sn * dV) + +# Watts +input_power = (3.0 / 2) * (Spe_tot + Spi_tot) * power_norm +input_pot = Sn_tot * Nnorm * rho_s ** 3 * 13.6 * 1.602e-19 * wci + +### Core input + +print("\n==== CORE ====") +print(f"Simulation time: {t_arr[-1]/wci * 1e3} ms") + +print("\n==== CORE ====") +print( + f"Input power : {input_power*1e-6:.2f} MW (heat) + {input_pot*1e-6:.2f} MW (potential energy)" +) +print(f"Peak input power density : {np.max(Spe + Spi)*pdens_norm * 1e-6:.2f} MW/m^3") +print(f"Electron heating: {(3./2)*(Spe_tot) * power_norm * 1e-6:.2f} MW") +print(f"Peak electron input power density : {np.max(Spe)*pdens_norm * 1e-6:.2f} MW/m^3") + +print( + f"Core electron->ion transfer: {np.sum(Wi * dV * mask) * power_norm * 1e-6:.2f} MW" +) +print( + f"Max core electron->ion transfer: {np.max(Wi * mask) * pdens_norm * 1e-6:.2f} MW/m^3" +) +print( + f"Min core electron->ion transfer: {np.min(Wi * mask) * pdens_norm * 1e-6:.2f} MW/m^3" +) + +zrad_core = np.sum(Rzrad * dV * mask) * power_norm +print(f"Core radiated impurity power: {zrad_core * 1e-6:.2f} MW") +print(f"Peak core impurity power: {np.max(Rzrad * mask)*pdens_norm * 1e-6:.2f} MW/m^3") + +hyd_core = np.sum(Rp * dV * mask) * power_norm +print(f"Core radiated hydrogen power: {hyd_core * 1e-6:.2f} MW") +print(f"Peak core hydrogen power: {np.max(Rp * mask)*pdens_norm * 1e-6:.2f} MW/m^3") + +print("\n==== SOL ====") + +sol_mask = 1 - mask + +print( + f"SOL radiated impurity power: {np.sum(Rzrad * dV * sol_mask) * power_norm * 1e-6:.2f} MW" +) +print( + f"SOL peak impurity power: {np.max(Rzrad * sol_mask) * pdens_norm * 1e-6:.2f} MW/m^3" +) +print( + f"SOL radiated hydrogen power: {np.sum(Rp * dV * sol_mask) * power_norm * 1e-6:.2f} MW" +) +print( + f"SOL peak hydrogen power: {np.max(Rp * sol_mask) * pdens_norm * 1e-6:.2f} MW/m^3" +) + +print( + f"SOL transfer ions->neut: {np.sum(Qi * dV * sol_mask) * power_norm * 1e-6:.2f} MW" +) +print( + f"SOL max ion->neutral transfer power: {np.max(Qi * sol_mask) * pdens_norm * 1e-6:.2f} MW/m^3" +) +print( + f"SOL min ion->neutral transfer power: {np.min(Qi * sol_mask) * pdens_norm * 1e-6:.2f} MW/m^3" +)