Error Estimation for Plasma Power Deposition in Tokamak First Wall Designs

Accurate computation of power deposition is expected to be critical for the successful design of plasma-facing components (PFCs) in reactor-scale magnetic confinement devices. The work provides analysis and model computations to help explain the high levels of accuracy obtainable using relatively coarse meshing of PFCs to treat cases with sharp-edged shadowing of one PFC by another. It emerges that even small misalignments of the surface grid with respect to the shadow edge or “terminator” may greatly improve the accuracy to which the total power to the tile is calculated, hence a strategy for error control emerges which emphasizes meshing with different constraint rather than different size of triangulation.


I. INTRODUCTION
E XPERIENCE with magnetically confined hot plasmas indicates that a significant fraction of the power may escape through turbulence effects toward the walls of the confinement device.The associated leakage mechanism is indeed likely to be exploited to allow for Helium "ash" removal in a fusion reactor.Further, although it is expected that in steady-state reactor operation, power in the escaping particles may be diffused by interaction with a layer of largely neutral gas, most reactor designs envisage startup and ramp-down situations where the hot plasma contacts the plasma-facing components (PFCs) of the device walls, and scenarios where the neutral "insulation" may be transiently ineffective.In either case, experiment indicates that the power directly deposited on the PFCs by the particles may be concentrated in a relatively small region.The turbulent nature of the mechanism makes it difficult to be precise about the size to be expected in a reactor-scale device-early estimates indicated that the power density could be damagingly intense unless the PFCs were carefully designed, and it remains a potential constraint on reactor operation.The author is with U.K. Atomic Energy Authority, Culham Science Centre, OX14 3DB Abingdon, U.K. (e-mail: wayne.arter@ieee.org).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TPS.2024.3399929.
Digital Object Identifier 10.1109/TPS.2024.3399929 Tools such as the SMARDDA-PFC program from the SMARDDA modules [1], [2], [3] are in routine use to support the design of PFCs in reactors based on the principle of the tokamak magnetic confinement device (see [4], [5], [6], [7]).Given the potentially critical nature of the code output, it is important, particularly when λ q , the width of the region of which power is lost, becomes narrow, to be sure that power is accurately deposited.Leaving aside questions of the accuracy of the physical model, for cases where λ q is of order 5-10 mm, it is a concern that the PFCs which are separately of meter-scale, are sufficiently accurately represented in the geometrical model.Customarily, the PFCs are represented as by a triangulation of their surfaces, which the preceding estimates suggest might each need millimeterscale sampling and thus contain a million triangles.For each triangle, a fieldline is traced using an adaptive Runge-Kutta integration, to determine whether it is shadowed by other triangles.Thus, even though the SMARDDA algorithm for testing for intersection is relatively sophisticated, cost scales at least as strong as the product of the number of triangles tested and the average cost of a fieldline integration.Computation times even for one PFC on a desktop machines may therefore easily exceed the 10 s or so regarded as a comfortable limit for an interactive design capability.Furthermore, most design questions normally involve a significant fraction of a first wall composed of 100-1000 PFCs.
This potentially overwhelming cost was anticipated by Arter et al. [2] in that a meshing interface was implemented that allowed for refinement of each triangle by up to two successive factors of four, wherein each inserted point was located "on" the geometry produced by the computer-aided design (CAD) system.It was planned that eventually refinement might proceed only in regions where the surface power density distribution Q varied rapidly.However, this was never implemented and the interface became neglected in favor of meshers that produced a single level of triangulation.The switch was possible because experience as illustrated in Section II indicated that, at least the total power P T deposited on the PFCs may be accurately accounted for on triangulations which are centimeter-scale.The accuracy to which the total power deposited on a PFC tile may be calculated by coarse meshes is upon consideration remarkable, when it is remembered that in case of a sharp-edged shadow thrown by one PFC on another, quadrature involves a discontinuous function.Under such circumstances, results are expected to be no better than first-order accurate (see Section III), thus for meshes with triangle side of order, the profile scalelength, order unity errors might naively be expected.
Since the present computational strategy does not allow for selective mesh refinement, computing P T is analogous to pixel-based computer graphics rendering problems, such as illustrated in Appendix A. However, in the PFC context, the "pixels" are triangular not rectangular and may have different sizes, shapes, and orientations.There is a good deal known about locating edges in pixellated images, edge information which could be exploited to produce more accurate quadrature, except that the literature deals with regular planar tessellations, and therefore does not seem easily transferable to more irregularly meshed 3-D surfaces.Moreover, the question of quadrature over pixels does not really arise in computer graphics, cf. the treatment of sampling in [8,Sec. 7].
The most relevant work found is that of Kamgar-Parsi and Kamgar-Parsi [9], who consider statistical measures of approximation accuracy only for differentiable functions, restricting attention to information from pixels forming regular planar tessellations with triangles, squares, and hexagons.Kamgar-Parsi and Kamgar-Parsi [9, Fig. 2] imply that triangles have poorer approximation properties than squares that in turn perform worse than hexagons.The current meshing interfaces for complicated surfaces remove any shape other than triangles from consideration, and Kamgar-Parsis' work indicates that it could be worthwhile to consider using them in a hexagonal pattern in addition to the Union Jack pattern previously employed [2] (see Section IV-A).
For these two tessellations with specially aligned shadow edges, Section IV-B describes the theoretical results that bear out the pessimistic estimates and serve a basis for cross-checking calculations made by the specially written code of Section V-A.The new IPROG software is capable of dealing with arbitrary configurations of terminator and triangulation and demonstrates that generically these dramatically reduce the error in P T .The reasons for the reduction is traced in Section V to a quasi-Monte-Carlo (QMC) sampling of Q effected by nonalignment.The concluding Section VI contains a summary.

II. PREVIOUS PFC COMPUTATIONS
Tests of the ability of SMARDDA-PFC to conserve total energy, equivalently total power in the model context, have been performed separately for many of the different tokamaks considered.The most detailed sensitivity studies were performed for the JET tokamak, which if not fully reactor-scale has a long history of experiment, hence the focus here on this device.A study denoted as 271.1 is representative and is described in detail in Appendix A. Although the Eich profile is used in the SMARDDA-PFC calculations, the profile parameters are such that to a fair approximation, the power deposition outside the last closed flux surface (LCFS) at midplane is an exponential with a fall-off length λ q = 17 mm.
The total power to be deposited is 10.2 MW, but to limit computational expense, attention was restricted to one 15 • segment of the divertor, whence 425 kW might have been expected.However, the smallest gap at midplane between the LCFS and the first wall is approximately 5 cm which is only 3λ q ; thus, a fraction e −3 ≈ 0.05 of the power is unavailable for assignment to fieldlines in the SMARDDA-PFC computation and a value of P T ≈ 400 kW is found in consequence (this truncation of the exponential profile would have to be treated by a linear scaling of Q values by ≈1.06 in the event that surface temperature rises were to be calculated but for present purposes, it is adequate simply to regard P T as being reduced).Plots of the power deposited on the divertor are of interest (see Fig. 1), where it will be seen from T6 in Fig. 1(b) that shadowing leads to sharp cutoffs both in the toroidal direction (decreasing y in the figure) and the radial direction corresponding to increasing x.Although the former is a commoner situation, the radial cutoff is also potentially relevant to reactor operation, and the more challenging because Q is rising exponentially fast at cutoff, hence this latter situation will be the focus of investigation herein.
For the specified equilibrium, the tile row handling peak power is that denoted T6, where triangle side-length h is estimated using ParaView [10] by sampling a number N of triangles on the exposed surface, calculating their total area A and using h = (A /(2N )) 1/2 .Four calculations were performed, and avoiding tabulation, the results expressed as pairs "(h in mm, P T in kW)" were (16 432), (10 390), (8406) and (3397).It will be seen that although h ≈ λ q , the SMARDDA-PFC-ParaView combination gives a value of P T within 10% of the likely true value, whereas naive estimates h /λ q of error in Section III are an order of magnitude larger.This unexpected accuracy is a feature of many SMARDDA-PFC calculations, and its explanation is the main aim of this work.

III. NAIVE ERROR ESTIMATES
SMARDDA-PFC's principal role is to decide whether a given triangle is shadowed by another one, for which the criterion is that the barycentre lies in the "lit" rather than the shadowed region.The code then assigns a single value of power density Q = Q i to each illuminated triangle i, from which the total power on a tile is calculated as a simple area-weighted sum of where A i the triangle area.The overall dimensions of each tile are of order 200 mm or more, meshed with as few as 200 triangles in the JET examples of Section II.ITER tiles have typical dimensions of order a meter, and this size was used in the computations of Section V-A as more reactorrelevant.
The discussion of Section II covers mesh sizes up to 16 mm, whereas the midplane Q-profile has an Eich profile, which falls off as R > R m increases at a rate close to exponential over a lengthscale λ q = 17 mm, and drops to zero in R < R m with a lengthscale of 1.1 mm.Both these lengths are, unless the discharge is limited, appreciably shorter than those which are realized on the PFCs, for two reasons, flux expansion and due to the effect of projection on the surface.
The result follows from elementary geometrical considerations given the power deposition formula, which say in the case of the exponential profile is given at midplane as where suffix "m" denotes a quantity evaluated at the largest major radius (R = R m ) of the LCFS.Elsewhere for where ψ is the local value of poloidal magnetic flux and B p = |∇ψ|/R is the strength of the poloidal magnetic field in Wb/rad.Suppose a surface with unit normal n to have an arc length coordinate s 2 in a vertical plane through the device major axis, i.e., a poloidal plane, so that where e φ is the unit vector in the toroidal direction, hence and thus where angle θ p is defined in the poloidal plane (see Fig. 2).Substituting in (3) for ψ − ψ m = ψ from ( 6), leads to a formula for the profile spreading length over the surface where it is important to note from Fig. 2 that the size of θ p may be much greater than θ , the angle between the total magnetic field and the surface normal.The product of the last two terms in (7) will be seen from the definition of B p to give a quotient of values of ψ at midplane and at the PFC surface, i.e., the flux expansion.In conventional tokamaks, this is often not much more than one, except near the X-points, as may be seen from use of the circular cylindrical formula for poloidal field where I φ is plasma current, and if (r, θ ) are cylindrical coordinates about the axis through the plasma center, R = R 0 + r cos θ so that so that the product R B p generally changes more slowly even than B p as a function of r .The projection effect may be seen from the flux plots to be relatively small for most of the PFCs where θ p ≈ 30 • translates into a factor of only 1/ sin(30 The preceding arguments suggest that λ f may nonetheless be of order triangle sidelength h = 16 mm, whence the simple estimate of first order accuracy∝ h /λ f still points to gross inaccuracy in the calculation of P T .The immediately following section serves to quantify and confirm this high level of error.

IV. MESH-ALIGNED SHADOWING
A. Patterns Fig. 3 shows the two triangular patterns of pixels used in the present study.Note that the "hexagon" pattern has in fact been generated by the shearing of a quadrilateral unit cell.There are well-known to be, in group theoretic terms, 17 wallpaper patterns (see [11,Sec. 26]), in which sense these two patterns are distinct, as the Union Jack is classified as p4mm whereas the hexagons area in class p3m1.The remaining 15 patterns are ruled out, either because they involve quadrilaterals, or because the relevant meshing capability is unavailable.Note that most tiles are not precisely rectangular because they form part of an axisymmetric device surface, hence even the "Union Jack" meshing of JET PFCs varies significantly from the idealization shown in Fig. 3, while the more usual Delaunay meshing may be very different from the regular hexagonal pattern drawn.

B. Theory
For the regular pixel patterns of Fig. 3 with the shadow edge aligned with the vertical, the integrated power deposited by the exponential profile may be calculated analytically, using the properties of a geometric progression (GP).It is convenient to work in terms of the horizontal coordinate normalized by λ f , i.e., x/λ f → x.Let the shadow edge be at x = d.Imagine for convenience and without loss of generality that |d| is close to zero.Each triangle has the same area so the quadrature equation ( 1) reduces to a summation.The horizontal dimension of the tile may be assumed sufficiently large that the power deposition Q is negligible at large x, and for definiteness suppose that the "height" or lateral y dimension of the tile is b.Analytic evaluation of the total power deposition per cell height gives simply 1) Union Jack Mesh: For the Union Jack mesh, each triangle has the same area A 1 = h 2 /8 = h 2 /2, but the triangle barycentres are not uniformly distributed in x, lying at |x| = h /3 = h/6 and |x| = 2h /3 = h/3, but not at distance |x| = h = h/2 from the left-hand edge of the pattern in Fig. 3(a).Periodicity in x implies that it may be assumed that −h/6 < d < h/3.Over a row of Union Jack flags, the quadrature sum I UJ (d) thus takes two different forms, depending whether the shadow boundary lies closer to edge or the center of the flag, precisely where h 3 and h 6 have been introduced to simplify expressions which include an offset d measured with respect to the left-hand edge in Fig. 3(a).Remembering the normalization of the x-coordinate by λ f , Q(x) = exp (−x) for the exponential Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
distribution.Each sum is then seen to consist of two GPs, thus which on rearranging give whence, for a sufficiently wide tile the GP sums give formulas Note that plotted as a function of d, two discontinuities in the quadrature I UJ (d) per unit cell.For small h, the integrals per unit height are approximately which simplify further as from which it follows that estimates of I range from 1−h/6 to 1 + h/6, to be compared with the analytic value I a = 1.Thus, although the relative error is, as expected O(h), the coefficient is favorable, so that inaccuracy, remembering the scaling of coordinates above, may be estimated as of order 1/8(h /λ f ) in the mean where h is the dimensional value of triangle side.More precisely, define for I (d) the mean absolute error as for a generic unit cell with twofold symmetry.For the Union Jack pattern, Hence, ϵ maUJ = 5h /36 ≈ 0.139 h .
2) Hexagonal Mesh: For the hexagonal mesh, each triangle has the same area A 2 = h 2 √ 3/16 = h 2 √ 3/4, and the triangle barycentres are now uniformly distributed in x, lying at x = 0, x = h /2 = h/4, x = h = h/2 and at x = 3h /2 = 3h/4.It may be assumed that 0 < d < h/4.Over a row of hexagons, the quadrature sum I HX (d) thus takes only the single form Remembering the normalization of the x-coordinate by λ f , Q(x) = exp −x for the exponential distribution, the sum over a single GP, thus whence, for a sufficiently wide tile, the GP sum formula yields per unit height Note that plotted as a function of d, a discontinuity every h/4 in the quadrature I HX (d).For small h, the integral per unit height is approximately e d−h/4 (1 + h/8), and since the minimum value of e d−h/4 ≈ 1 − h/4, the range of values is approximately 1 − h/8 < I < 1 + h/8 to be compared with the analytic value I a = 1.Thus, although the relative error is, as expected O(h), the coefficient is favorable, so that inaccuracy, remembering the scaling of coordinates above, may be estimated as of order 1/8(h /λ f ) in the mean where h is the dimensional value of triangle side.Indeed, using the definition equation (23), since there is fourfold symmetry Hence, ϵ maHX = h /8 = 0.125 h .

C. Quadrature by Trapezoidal Rule
It will be seen that in the presence of a boundary located between sample points, it is better to use the simple sum over central points equation (1) rather the trapezoidal rule for integrating a function.Consider first the simple sum applied to the problem of equidistant sample points separated by a distance h, then the first point within the illuminated region is subjected to a power density value of ), so the approximation to the integral is whence, for a sufficiently wide tile, the GP sum formula yields per unit height Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
and on averaging over −h/2 < d < h/2, the simple sum I SP is unity to order h 2 .Contrast this result with the fact that the trapezoidal rule simply reduces the (leading) factor of unity in (29) to one-half, supposing that it is most consistent simply to omit the contribution of any h-interval in which the shadow edge lies.Hence, the trapezoidal rule in the mean will give a smaller value than the expected result by a factor of 1 − (h/2)e −h/2 .

V. ARBITRARY SHADOWING A. Computations
The repository SMARDDA-QPROG [12], designed primarily to provide a test harness for new SMARDDA modules, contains a suitable pair of skeleton objects for producing code SMARDDA-IPROG to test the quadrature or "Integration" properties of different mesh arrangements.In the present case, although it is unlikely that the new objects would be incorporated in the core SMARDDA software, it is still helpful to work within software layouts and conventions familiar to SMARDDA developers.
For the current work, following edits to qprog.txt and bigobj.txt, the set-up script was invoked by ./qprog.bash -l QPROG=sumtot Q=i \ STR="integral_estimates" \ BIGOBJ=sum BO=su BSTR="box_and_function." The script edits qprog.f90 to give a new program sumtot.f90 with control module icontrol, and edits bigobj_m.f90to generate an associated module sumtot_m.f90.Controls of type inumerics are invoked in icontrol, using inputs described in qprog.txt.The script further edits bigobj_m.f90to generate a module object sum_m.f90from the description provided in bigobj.txt,plus a second type (sunumerics) to transfer parameter inputs described in qprog.txt.
Thereafter, there were user edits to sumtot_m.f90,formerly bigobj_m.f90,as follows: 1) sumtot_dia to call object diagnostics control routine sum_dia; 2) sumtot_solve to select pixel pattern and loop over mesh displacements; 3) sumtot_write to list key input parameters to the.out file; 4) sumtot_writeg to call gnuplot output control routine sum_writeg.And user edits to sum_m.f90,formerly bigobj_m.f90,as follows.
1) Add to sum_solve to implement the quadratures.
2) Duplicate sum_solve as sum_solvehex to implement the quadratures over hexagons.3) Add to sum_dia for new output to.log file.4) Add to sum_writeg and copy to sum_writehexg to output picture of tessellation to.gnu file.Consistent with SMARDDA practice, SI units were used for lengths, while angles were assumed to be in degrees (usually radians could also be specified subject to use of an input switch).The "unit cell" size h was used in input, rather than triangle sidelength.The exponential profile, with its discontinuity at the terminator, was used throughout.Each tile, consistent with ITER dimensions, was assumed to be a × b = 1 × 1 m, with λ f = 0.1 m, hence there were approximately ten e-folding lengths across the tile in x, so that the total power deposited would be expected to deviate negligibly from unity (note of course that values O(1) MW, would be expected in practice, but the total power may be scaled without affecting relative accuracy of computation).The code works by assuming that the shadow boundary is fixed to lie along the y-axis, and the triangle pattern is rotated with respect to y by an angle β (see Fig. 4).For perfect alignment of the shadow boundary with the y-axis (i.e., normal to one cell boundary), lateral dimension b is irrelevant, but when there is misalignment by an angle β > 0, then the number N y of cells in y is important, N y = b cos β/ h = 2b cos β/ h for the Union Jack pattern, and N y = 2b cos β/( √ 3h) for hexagons.Although this implies that there are approximately 14% more hexagon rows, this is not felt likely to be significant in the present context.
The variation of the computed integral I as d varies for β = 0 is plotted in Fig. 5 for the Union Jack meshing and in Fig. 6 for hexagons, showing both qualitative and quantitative agreement with the analysis presented in Section IV-B.From the data produced by these computations are derived the maximum relative error as h and β vary for the most significant calculations performed using SMARDDA-IPROG, which are listed in Table I.For the Union Jack tessellation, the errors are visualized as a function of d for different mesh spacings in Fig. 7 and for different mesh orientations in Fig. 8.The results at different orientations of the hexagonal mesh are shown in Fig. 9.The last is notable in that the rotation through 30 • aligns Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the shadow boundary with nonhorizontal hexagon edges, and recreates the nonuniform sampling in x found for the Union Jack case.
The small size of the error for many of the tilted cases is understood because there are now up to N b different sample values corresponding to each of the points in the "vertical" cases.Provided the tilt is such that a region h distant from the terminator is well sampled, i.e., b| sin β| ≥ h , equivalently β ≥ h /b (but not a special value such as β = 30 • in  the hexagonal mesh), the resulting sampling of Q(x) will generically be QMC, analogous to lattice rule evaluation of the integral [13] (see Appendix A).For such rules, an error scaling as (ln N b ) 2 /N b is expected, and since N b varies from N b = 32 when there are three points in λ f , to N b = 200 when there are 20, reductions of up to two orders of magnitude in the quadrature error found at β = 0 are expected and confirmed by Table I.The minimum value of β needed for the increased accuracy is seen to be modest, as 1/32 rad corresponds to 2 • .Of course, for small β, there are also "magic" values of the tilt which effectively provide uniform sampling at a separation of N t h sin β, where N t is a small integer, giving an error That quasi-random sampling is effected by the tilt, is supported by 1-D calculations comparing both Monte-Carlo and QMC evaluation of the integral of the decaying exponential.The random numbers for the Monte-Carlo sampling were provided by the "Luxury" pseudorandom number generator of Marsaglia and Zaman, implemented by James [14].Quasi-random numbers were provided as a Halton sequence (see [15]).The source code for both generators was patched into SMARDDA-IPROG to provide a comprehensive record.For the bulk of the tilted calculations, N b = 5000, hence this number of samples was used in the numerical experiments.For both the pseudorandom and the quasi-random generator, sequences of ten estimates were made of the integrated power and its average deviation, minimum deviation and maximum deviation from unity were recorded.For the average deviation, the pseudorandom calculations gave a value of 0.0130 compared to 0.0007 for quasirandom, and the maximum deviations were, respectively, 0.04 and 0.003.Comparison with the values in Table I indicates that the quasirandom statistics give much the better fit (it is worth noting that the estimate of expected QMC error in Appendix A appears to be particularly unfavorable relative to MC error, since the theory predicts only a 70% reduction, not the factor of approximately ten deducible from Table I).

VI. CONCLUSION
The work has demonstrated that relatively coarse meshes, using triangles of side h ≈ λ f /3 may give tile-integrated powers P T accurate to within 1% provided that the following holds.
1) The mesh is not specially aligned with the shadow boundary, so that edges deviate by at least an amount of order h .2) The shadow boundary is of order ten or more h in length.Overall, there appears to be little difference in accuracy between the Union Jack and the other triangular mesh pattern.However, it is evident that a good strategy for verifying convergence is to produce meshes with different constraints, most obviously to drive mesh alignments with differently oriented vector fields.In practice, relatively small changes in the initial node spacing along an edge or in mesh quality criteria such as the minimum allowed angle, will often be seen to produce suitable effects in meshes on 3-D PFC surfaces.Shadow edges which curve differently to the pattern of surface tessellation are also similarly expected to be beneficial.
The location and maximum value of Q is less important for the reason that it is expected that it is maximum value of temperature T max is the most important quantity as affecting the material properties of the surface.Computed T max is sensitive to the properties of the (probably) metallic material used in its construction, so that proper consideration of local accuracy requires study of a coupled problem which is outside the scope of the present work.Arguably for the relatively specialized case of a strike-point (effectively a shadow edge) sweeping over 3 cm radially, conserving total Q is more critical than locating T max precisely.There is any event, the practical constraint that thermal calculations normally involve a 3-D mesh which computational costs will demand have edge lengths much greater than a 2-D triangulation.Nonetheless, elementary considerations imply, at least in the case of sharpedged shadows, that similar conditions as for P T apply to achieve improved accuracy for the maximum Q computed value on a tile.The extremum value on the mesh will be given by the triangle with its center nearest the shadow edge, so it is important to have a good spread of central point distance about the terminator.
It is necessary to conclude with the warning that not all shadowing is as simple as has been assumed, and particularly for one-off calculations, a successive refinement strategy should always be considered.One common example where small triangles may be essential to capture key features of the deposition is a case where a PFC tile face just fails to shadow, thereby slightly exposing, the edge of an adjacent tile.More subtle effects are conceivable.

APPENDIX A PROVENANCING OF JET WORK
The presence in the Appendix rather than the main text, of details of provenance, reflects the fact that the detailed reproduction of a particular JET case is not crucial in the present context where a generic effect is claimed.The main provenancing of the JET studies mentioned is through a git repository local to UKAEA, a clone of which can be provided on request.The major subdirectories of the repository all have README files to describe their content and in some cases, point to READMEs in their subdirectories.To predict Q deposited on the first wall, the SMARDDA-PFC software requires three main pieces of information, namely: 1) magnetic equilibrium; 2) parameters describing power deposition; 3) geometry of first wall.The magnetic equilibria used in the JET studies are saved as eqdsk files in subdirectory Data/Equilibrium of the repo and the file provenance.txtdescribes how they were obtained.File g_p90271_t49.0_129× 129.eqdsk defines the equilibrium 271.1, which has toroidal magnetic field B T = 2.84 T and I ≈ 2.4 MA.
The parameters λ q = 0.017 m S = 0.0011 m and P tot = 10.2MW describing power deposition for the JET pulse #90271 at T = 49 s were provided by Silburn et al. [16] of UKAEA using software written by him for parametric fitting to JET thermal data.The supply of the JET divertor geometry is described in subdirectory Data/STP.
For completeness, it is recorded that JET pulse #90271 has an input neutral beam power of 20 MW, and is a strike-point sweeping experiment, i.e., the current varies by a percent or so such that the flux contours move back and forth by approximately 3 cm with a frequency of 4 Hz.Additionally also, it is noted that the Eich power deposition profile employed has the form at midplane, in the above parameters, of (31) Fig. 10.Related computer graphics sampling problem.Each square of the lattice corresponds to a pixel used to render the image, and in the simplest algorithm, the pixel will be regarded as "full" of circle if its center lies within the circle boundary.
where s = (S/2λ q ), and x = ( ψ/R m B pm ), and the erfc function is defined as It is mathematically striaghtforward to normalize the profile so that P tot leaves the midplane, and to show that the profile rises from zero over a distance of order S to its peak and thereafter falls off approximately exponentially at a rate λ q if S ≪ λ q .
APPENDIX B RANDOM AND QUASI-RANDOM NUMBERS Monte-Carlo methods have a long and excellent history in many areas of science and technology, including plasma physics.Nevertheless, Monte-Carlo methods have the drawback that their error ϵ N is usually where N r is the number of randomized samples used in the calculation.Hence, compared to a deterministic method of order p accuracy in d dimensions, with they may exhibit very much slower when p > d/2.Conversely, this is the reason why Monte-Carlo is preferred for integration of functions of many variables, i.e., high d-dimensionality integrals.
There are other situations where order p methods may be of less utility, notably where the quantity being calculated is badly behaved.An example is integration when either the integrand or the domain of integration is not smooth, so that theoretical degree of convergence is not achievable, instead where p ′ < p.When estimating area, typically p ′ = 1, as may be seen from the problem of estimating the area of a circle using a uniform square lattice of integers with side 2n L , so N r = (2n L ) 2 .Suppose that the circle is comparable in size to the entire lattice, but fitting within it, i.e., with radius a 0 = c 1 n L , where c 1 < 1 (see Fig. 10).Assuming a 0 ≫ 1, the circumference of the circle intersects of order 2πa 0 lattice cells.The average error associated with each boundary cell may be estimated as 0.5, hence the total error in area as πa 0 , and the relative error in area is 1/a 0 .Since N r = (2a 0 /c 1 ) 2 , it follows that ϵ N ∝ 1/N r 1/2 in agreement with (35).The above difficulties with both Monte-Carlo and deterministic methods have motivated the development of QMC methods [13], which have typically the property The principal source for QMC methods is the book by Niederreiter [13], which however is highly mathematical in tone.
A good way to understanding is to start by considering the subject of quadrature, specifically looking at an integral over two dimensions representing the area of a peculiarly shaped finite object.Suppose the object fits in the unit square.Mathematically, it is defined by 1, x within and on the object 0, otherwise.
The required area integral (and indeed for general integrands) is then where dτ is a volume element in d-dimensional space (so dτ = d xdy in 2-D).The Monte-Carlo approximation to the above integral equation (37) is defined as where x n are vectors of random variables on the unit interval [0, 1).Any textbook (e.g., [13, Sec.1]) will explain that the error in the integral is where σ ( f ) is the standard deviation of f in the region of integration.The error in the approximation equation (38) may be expressed in terms of the "point set discrepancy" or discrepancy for short (the discrepancy for randomly distributed points is ∝ 1/ √ N r , hence (39)).Unlike deterministic integration methods where the quadrature points have fixed locations and so the error depends only on the form of the function f , for Monte-Carlo the error depends on how the sample points are distributed.Hence, the use of discrepancy to measure error, since discrepancy is measured with respect to a certain set of subintervals, specifically 2-D rectangles within the unit square.Precisely, consider the errors made in calculating the area of the object within each rectangle in the set, then discrepancy measures the largest such error.There are in fact two closely related measures of discrepancy, depending on whether the subintervals all include the origin ("star discrepancy") or whether the set consists of all possible rectangles ("extreme discrepancy").The distinction is not of great importance because the two measures lie within a factor of two of each other.
In 1-D, it is easy to show that the discrepancy, however, measured is least D N ∝ 1/N r when the sample points x n are uniformly distributed.The key fact is that there exist sets of points ("low discrepancy sequences") for which the discrepancy does not increase very much as the number of dimensions increases.
The simplest of these sets to describe is that due to Halton.In 1-D, it is identical to a van der Corput sequence, which involves generating numbers on the unit interval, using the reversed bit patterns of the positive integers.It is best illustrated by example.Thus, 2 has the binary representation 10, so the second element in the van der Corput sequence is .01 2 or 1/4, 3 = 11 2 , so the third element in the van der Corput sequence is .11 2 or 3/4, 4 = 100 2 , so the fourth element in the van der Corput sequence is .001 2 or 1/8.Hence, the first seven elements of the van der Corput sequence are (in eighths) 4, 2, 6, 1, 5, 3, 7. (40) It will be seen that there is a sort of fractal pattern about the above distribution.Van der Corput sequences may be defined for any prime b, by representing the integers in the base b, then using the reversed representation as above to generate values on the unit interval.The Halton sequence in 2-D contains pairs of numbers, the first in the nth pair being the nth element from a van der Corput sequence with base 2 and the second being the corresponding element in a base 3 van der Corput sequence.In fact, any two distinct primes could be used, and the generalization to many dimensions should be obvious.The discrepancy of the Halton sequence in 2-D is bounded by a formula which may be approximated as where A 2 = 0.66.It will be seen that this is smaller than the Monte-Carlo value of ϵ N = 1/ √ N r for N r > 1000.It is therefore also competitive with uniform sampling on a rectangular lattice of N r 2-D points, which in general gives an error proportional to lattice spacing, i.e., 1/ √ N r .

Manuscript received 9
January 2024; revised 9 April 2024; accepted 8 May 2024.Date of publication 20 May 2024; date of current version 14 August 2024.This work was supported in part by the European Union via the Euratom Research and Training Programme under the framework of the EUROfusion Consortium under Grant 101052200 (EUROfusion) and in part by the Engineering and Physical Sciences Research Council (EPSRC) under Grant EP/W006839/1.The review of this article was arranged by Senior Editor S. J. Gitomer.

Fig. 1 .
Fig. 1.Power deposition Q in MWm −2 for a 15 • toroidal segment of JET divertor in magnetic field eqid = 271.1 as described in the text, for (a) coarsest and (b) finest griddings considered.The tile rows are labeled in (a).The plasma is sited above the divertor geometry as indicated by Fig. 2.

Fig. 2 .
Fig. 2. Flux surfaces for eqid = 271.1,illustrating the angle θ p .This is the angle at intersection between made between the outline or silhouette of the JET PFCs shown in green and the purple contours of flux as proxies for projections of fieldlines onto a poloidal plane.Hot plasma lies in the region above and bounded by the lowest concave upward flux contour, which represents the LCFS ψ = ψ m , passing through the midplane Z = 0 at R = R m .Typical values of θ p are seen to be θ p ≃ 30 • , ranging up to approx.60 • at bottom left.

Fig. 3 .
Fig. 3. Illustration of pixel patterns.(a) "Union Jack" pattern.(b) "Hexagonal" pattern.Dots correspond to the barycentres of the triangles.Note in both cases that the "unit cell" of the pattern has a horizontal dimension h = 2h , where h is the length of a side of the triangle.Coordinate x increases from x = 0 at far left, and may be identified with the direction of flux variation, implying the vertical coordinate corresponds to the toroidal direction.

Fig. 4 .
Fig. 4. Illustration of a misaligned shadow on "hexagonal" pattern.The horizontal x-direction corresponds to the radial coordinate directed approximately in the X -direction of Fig. 1, and the vertical to the toroidal, approximately the Y -direction of Fig. 1.

Fig. 5 .
Fig. 5. Integral as a function of mesh displacement d for the mesh-aligned shadow boundary.Union Jack pixellation, h = 0.01.

Fig. 6 .
Fig. 6.Integral as a function of mesh displacement d for the mesh-aligned shadow boundary.Hexagonal tessellation, h = 0.005.

Fig. 7 .
Fig. 7. Integral as a function of mesh displacement d for different values of mesh resolution listed on the figure, β = 2.5 • .Union Jack pixellation.

Fig. 8 .
Fig. 8. Integral as a function of mesh displacement d for different values of mesh rotation listed on the figure.Union Jack pixellation, h = 0.03.

Fig. 9 .
Fig. 9. Integral as a function of mesh displacement d for different values of mesh rotation listed on the figure.Hexagonal tessellation, h = 0.005.