NTCC module -- kdsaw -- Kadomtsev Sawtooth Model with Porcelli Upgrade. ======================================================================= D. McCune (dmc) 29 Jan 2002 Update: dmc 15 Nov 2005 -- Controls added for Porcelli "partial reconnection" mixing model; documentation updated. Abstract: This module implements the Kadomtsev Sawtooth Mixing Model[1] and the Porcelli Sawtooth Mixing Model[2] for axisymmetric tokamaks in general geometry. The module is initialized by calls to "kdsaw_pre_init" and to "kdsaw_init", the former specifying the mixing model selection, and the latter providing a grid, the iota(bar) (1/q) profile, and zone volumes; kdsaw_init computes the Kadomtsev or Porcelli mapping and stores it in the module's internal memory. Subsequent calls are used to extract specific information about the sawtooth mixing. The grid passed to kdsaw_init is specified in normalized square-root-toroidal-flux, a coordinate which functions in a way strictly analogous to the minor radial coordinate "r" in Kadomtsev's original paper, but is suitable for general geometry. The grid is given explicitly and need not be evenly spaced. The current profile is expressed as (1/q), and the zone volumes (in any suitable units) are used to adjust the mapping such that a change in a particle or energy density profile according to the ideal Kadomtsev mapping volume integrates to precisely zero: typically, a peaked density profile is flattened, with the density reduced in the center and raised at the edge of the mixing region. A user-adjustable current profile remapping is also provided. Both REAL and REAL*8 interfaces are provided to all routines in the module. [1] Sov. J. Plasma Phys. Vol. 1 No. 5 Sept. Oct. 1975, p. 389. [2] Porcelli et al., Plasma Physics & Controlled Fusion 38 (1996) p. 2163, (See section 4). ------------------------------------------------ Note: the Porcelli sawtooth trigger model is NOT included in this software. A separate NTCC module is expected to be supplied for the sawtooth trigger calculation. KDSAW calculates mixing only. ------------------------------------------------ Kadomtsev mixing vs. Porcelli mixing: the details are provided in the references, but summary of the mixing characteristics are given here. The coordinate x refers to normalized sqrt(toroidal flux); x=0 means the magnetic axis, x=1 would refer to the plasma boundary. For both models, there is no "instantaneous" (i.e. MHD timescale) change in the plasma current profile or plasma parameters except inside the sawtooth mixing region. Let x_q1 be the x location of the q=1 surface which for purposes of this discussion is assumed to be unique; let x_mix > x_q1 be the outer limit of the sawtooth mixing region. ---------- Kadomtsev current mixing: within the mixing region, the q profile value is increased in the region 0 < x < x_q1, and reduced in the region x_q1 < x < x_mix, and somewhat flattened. The post sawtooth q profile has q(0)=1 and q'(x)>0, q(x)>1, for x>0, out to x_mix, where a "current sheet" (discontinuity in the q profile) forms, since the pre- and post-sawtooth q profile values are assumed unchanged outside the mixing region. Some experimental evidence suggests q<1 persisting after sawtooth crashes; this was part of the motivation for the later development of the Porcelli model. ---------- Kadomtsev plasma mixing: in the standard Kadomtsev model, the plasma at the pre-sawtooth q=1 surface is transported to the axis x=0. The pre-sawtooth q=1 surface represents the maximum value of a pre-sawtooth helical flux function Psi_1(x) = [poloidal magnetic flux from 0 to x] - [poloidal magnetic flux from 0 to x if q=1 everywhere]. Away from the pre-sawtooth q=1 surface, the plasma is zoned by matching value ranges of Psi_1, giving two plasma regions one inside and one outside the q=1 surface for each pair of Psi_1 values. The plasma from these two regions is mixed and transported to a region subtending equal toroidal flux, with zone pairs corresponding to higher pre-sawtooth Psi_1 mapped nearer to the axis, and pairs with lower Psi_1 values mapped towards the edge of the mixing region. The net effect of this is generally to flatten the plasma temperature and density profiles but there can be residual structure depending on the details of the pre-sawtooth profiles. KDSAW allows the option of "ergodic" mixing instead, in which all plasma is redistributed evenly throughout the mixing region, resulting in totally flat post-sawtooth particle and energy density profiles. ---------- The Porcelli model depends on an additional parameter Fp, which is the "fractional q=1 island width" prior to the sawtooth crash. An inner boundary of the q=1 island is defined: x_min = (1-Fp)*x_q1 where Fp=1 implies x_min=0, which would result in mixing similar to the Kadomtsev model. In the Porcelli model, the pre-sawtooth Psi_1(x_min) is found; the outer boundary of the q=1 island and the mixing region, x_mix, is found by the condition Psi_1(x_mix)=Psi_1(x_min). In the Kadomtsev model, by contrast, Psi_1(x_mix)=Psi_1(0)=0. The Porcelli mixing region is smaller: the x_mix of the Porcelli model is less than the x_mix of Kadomtsev by an amount depending on the value of Fp. The Porcelli model considers the region x_min <= x <= x_mix to be the q=1 island which gets "reconnected" during the sawtooth crash. The inside axial region 0 <= x < x_min is not included in this reconnection, but is considered to go through a period of ergodic field lines and rapid plasma mixing nonetheless. Porcelli current mixing: for x_min <= x <= x_mix, q = 1 is imposed. for 0 <= x < x_min, q = q[pre-sawtooth](x_min) is imposed. So the q profile becomes a two-step function within the mixing region, with discontinuities or "current sheets" at x_min and x_mix. (Subsequent application of a poloidal field diffusion equation will result in rapid decay or spreading of the current sheets). Porcelli plasma mixing: The plasma is redistributed evenly but separately within the two subregions, resulting in two-step particle and energy density profiles after the sawtooth crash. KDSAW allows the option of "ergodic" mixing over the ENTIRE mixing region instead, in which all plasma is redistributed evenly throughout, resulting in flat post-sawtooth temperature and density profiles (with no step at x=x_min). This choice (ion_mix=4) does not affect current mixing. One of the four KDSAW mixing options are chosen by means of a call to KDSAW_PRE_INIT: ion_mix zf_reconnect 1 (ignored) Standard Kadomtsev mixing (default). 2 (ignored) Kadomtsev current mixing; temperature and density profiles flat post-sawtooth. 3 =Fp Standard Porcelli mixing. 4 =Fp Porcelli current mixing; temperature and density profiles flat post-sawtooth. ------------------------------------------------ Current sheet smoothing: The Kadomtsev mixing model produces 1 current sheet at x=x_mix; the Porcelli model produces 2 current sheets at x=x_min and x=x_mix. For finite difference numerical codes, the effective thickness of the sheet will be determined by the grid size. This means that the dynamical details of the response of poloidal field diffusion post-sawtooth will be affected somewhat by grid size. To reduce this numerical effect, it is possible to specify a finite width (in sqrt toroidal flux space) over which the current sheet(s) is (are) to be smoothed out. See the description of subroutine KDSAW_DXI_QJUMP below. The default (if KDSAW_DXI_QJUMP is never called) is for no such smoothing. ================================================ ================================================ Test program: The program "kdsaw_test" demonstrates the use of kdsaw subroutines. It requires no input file. To test the module, build the program and unix> kdsaw_test > test.output unix> diff test.output kdsaw_test.output Differences, reflecting slight variations of floating point arithmetic from one computing platform to another, should be minimal. The program "kdsaw_test" exercises (and demonstrates) all mixing options; current sheets are not smoothed. ================================================ ================================================ Subroutine call interface: A fortran-77 style interface is provided. The REAL*8 interface is described. Append "_r4" to get the REAL interface (e.g. use "kdsaw_init_r4" instead of "kdsaw_init" if REAL arguments are to be provided. Notes: (1) If anything other than pure Kadomtsev mixing is desired, execute "call kdsaw_pre_init" first, to set the Porcelli mixing controls. (2) Next, execute "call kdsaw_init"! All other calls refer to data that was input in a preceding kdsaw_init call. Each "kdsaw_init" call updates the Kadomtsev/Porcelli sawtooth map. ----------------------------- **** (optional) set Porcelli mixing controls **** Default values (used if KDSAW_PRE_INIT is never called) correspond to: zf_reconnect = 1.0 ! real*8 ion_mix = 1 once the defaults are modified, the new values stay in effect until the next kdsaw_pre_init call. However, it is recommended to *always* call kdsaw_pre_init prior to kdsaw_init, in codes that support non-Kadomtsev sawtooth mixing as an option. SUBROUTINE KDSAW_PRE_INIT(zf_reconnect,ion_mix) real*8, intent(in) :: zf_reconnect ! reconnection fraction (Fp) integer, intent(in) :: ion_mix ! plasma species mixing option The meanings are as described above in the section describing the Kadomtsev and Porcelli mixing options. KDSAW_PRE_INIT reports no errors or messages. The following value ranges are silently enforced with min & max function calls: 1 <= ion_mix <= 4 ..and.. 0.0 <= zf_reconnect <= 1.0 NOTE: if the Porcelli model is used, current sheets, i.e. discontinuous jumps in the q profile, will be produced. In practice the current sheet width will be spread out over a grid spacing in real codes. This might cause problems when comparing sensitivity of details of code dynamics to grid size; therefore, a non-zero current sheet width may be specified via a separate call to "kdsaw_dxi_qjump"... The routine "kdsaw_pre_init" must be called *prior* to "kdsaw_init" to have an effect. ----------------------------- **** (optional) set current sheet width **** Default value: dxi_qjump = 0.0 ! real*8 SUBROUTINE KDSAW_DXI_QJUMP(dxi_qjump) real*8, intent(in) :: dxi_qjump ! current sheet width in xi space ! (i.e. sqrt toroidal flux normalized). Apply a smoothing operation of half-width dxi_qjump at each current sheet (i.e. the discontinuous changes in the post-sawtooth q profile at locations corresponding to the inner and outer boundaries of the q=1 magnetic island. For Kadomtsev mixing, only the outer boundary exists. The maximum allowed value if dxi_qjump is 0.1 -- this is silently enforced. Once the default is modified, the new value stays in effect until the next kdsaw_dxi_qjump call. However, it is recommended to *always* call kdsaw_dxi_qjump prior to kdsaw_init, in codes that support current sheet smoothing as an option. This routine must be called *prior* to "kdsaw_init" to have an effect. ----------------------------- ----------------------------- ----------------------------- **** compute Kadomtsev/Porcelli sawtooth map **** **** subroutine kdsaw_init **** SUBROUTINE KDSAW_INIT(ilun_msg,izones,zxi,ziotb,zdvol,iimix,ierr) ! input data: integer, intent(in) :: ilun_msg ! fortran i/o channel for error messages integer, intent(in) :: izones ! number of plasma zones real*8, dimension(izones+1) :: zxi ! zone *bdy* locations in normalized sqrt(tor.flux) ! in strict ascending order with... ! zxi(1)=0, zxi(izones+1)=1 expected ! zxi(j) is inner bdy, zxi(j+1) outer bdy of j'th zone real*8, dimension(izones+1) :: ziotb ! iota(bar) (1/q) profile at zone bdys ! These are pre-sawtooth values. real*8, dimension(izones) :: zdvol ! "volume" of the zones ! in arbitrary units, only the relative values ! amongst zones are important ! output data: integer, intent(out) :: iimix ! iimix gives index to last zone involved in ! Kadomtsev mixing, or 0 if none (not an error). integer, intent(out) :: ierr ! error flag-- error in input data, ! such as, izones.lt.5, or, zxi out of order. Note: zxi need not be evenly spaced, but it is crucial that its values correspond to sqrt(Phi/Philim) where Phi = enclosed toroidal flux. ----------------------------- **** retrieve mixing region information **** **** subroutine kdsaw_xmix **** SUBROUTINE KDSAW_XMIX(zxmin,zxq1,zxmix) real*8, intent(out) :: zxmin ! x_min -- inner q=1 island boundary real*8, intent(out) :: zxq1 ! x_q1 -- pre-sawtooth q=1 surface location real*8, intent(out) :: zxmix ! x_mix -- outer q=1 island boundary; outer ! limit of sawtooth mixing. If Porcelli mixing is used with Fp=0 (not recommended), zxmin=zxq1=zxmix; otherwise zxmin < zxq1 < zxmix is to be expected. For Kadomtsev options, zxmin=0.0 will be returned. The higher the value of Fp, the larger zxmix and the small zxmin will be (i.e. the larger the sawtooth mixing region). ----------------------------- **** compute post-sawtooth iota(bar) profile with ad hoc adjustment **** **** subroutine kdsaw_iotb **** SUBROUTINE KDSAW_IOTB(ziotb,zmix,ierr) ! return the post-sawtooth iota(bar) profile ! call this only after kdsaw_init has been called. ! user control: ! mixing parameter: zmix=1 ==> return fully mixed post-sawtooth profile ! zmix=0 ==> return pre-sawtooth profile unchanged ! 0 return ! (1-zmix)*[pre-sawtooth profile] + ! (zmix)*[post-sawtooth profile] ! ! CAUTION: zmix is NOT the control for Porcelli partial reconnection; it ! is a means for further reduction of the current mixing if desired. ! ! warning issued if zmix out of range; zmixl = max(0,min(1,zmix)) is used. ! ! note: dimension of ziotb array is equal to "izones+1" where "izones" ! was the dimensioning argument provided in an earlier call to kdsaw_init. real*8, intent(out) :: ziotb(izones+1) ! (mixed) post-sawtooth iota(bar) profile real*8, intent(in) :: zmix ! mixing parameter in range [0,1] integer, intent(out) :: ierr ! completion code, ierr=0 means OK ----------------------------- **** compute change in a density profile **** **** subroutine kdsaw_deln **** SUBROUTINE KDSAW_DELN(zpre,zdel,ierr) C C USING THE MIXING FACTORS CALCULATED IN THE SAWMIX ROUTINE, C EVALUATE THE EXPECTED CHANGE ZDEL IN PRE-SAWTOOTH DENSITY C PROFILE ZPRE DUE TO THE SAWTOOTH (ideal Kadomtsev/Porcelli mixing). C C THE VOLUME INTEGRAL ACROSS THE WHOLE PLASMA OF ZDEL WILL COME C OUT TO ZERO C C Note: the dimensioning of the arrays is determined by the argument C "izones" in a prior call to kdsaw_init. C real*8, intent(in) :: zpre(izones) ! input density profile real*8, intent(out) :: zdel(izones) ! ideal Kadomtsev change profile C integer, intent(out) :: ierr ! completion code, 0=OK ----------------------------- **** given a random number, choose post-sawtooth zone of particle starting in a particular pre-sawtooth zone. **** subroutine kdsaw_find **** SUBROUTINE KDSAW_FIND(zselect,izone_in,izone_out,zfrac) ! given an input zone and a (random) number between 0 and 1, select ! an output zone, according to the Kadomtsev redistribution mapping real*8,intent(in) :: zselect ! a number btw 0 and 1 integer,intent(in) :: izone_in ! pre-sawtooth originating zone integer, intent(out) :: izone_out ! post-sawtooth destination zone selected real*8,intent(out) :: zfrac ! location within post-sawtooth zone ! zfrac=0 corresponds to inner boundary; 1 corresponds to outer ! boundary of mapped-to zone. Note: if the input zone index indicates the last zone which has non-null intersection with the mixing region (i.e. a zone which is partly inside and partly outside the mixing region), the source point being mapped from is assumed to be inside the mixing region. For any call, if izone_out is the last zone in the mixing region, the maximum value of zfrac will reflect the fraction of the last zone that is in the mixing region. ----------------------------- **** given a post-sawtooth zone index show what fraction of material from which pre-sawtooth zones were mixed to form the post-sawtooth zone's contents. **** **** subroutine kdsaw_from **** SUBROUTINE KDSAW_FROM(izone,inzmax,inzact,ifrom,zfrom,ierr) ! give statistics on what fraction of material from which pre-sawtooth ! zones constitute the material in a selected post-sawtooth zone. integer, intent(in) :: izone ! selected post-sawtooth zone integer, intent(in) :: inzmax ! max number of pre-sawtooth zones integer, intent(out) :: inzact ! actual number of pre-sawtooth zones integer, intent(out) :: ifrom(inzmax) ! pre-sawtooth zone indices real*8, intent(out) :: zfrom(inzmax) ! pre-sawtooth material fractions integer, intent(out) :: ierr ! error code (0=OK) ----------------------------- **** given a pre-sawtooth zone index, show to which post-sawtooth zones the pre-sawtooth zone's material is distributed. **** subroutine kdsaw_to **** SUBROUTINE KDSAW_TO(izone,inzmax,inzact,ito,zto,ierr) ! give distribution of material of pre-sawtooth zone "izone" to ! affected post-sawtooth zones. integer, intent(in) :: izone ! selected pre-sawtooth zone integer, intent(in) :: inzmax ! max number of post-sawtooth zones integer, intent(out) :: inzact ! actual number of post-sawtooth zones integer, intent(out) :: ito(inzmax) ! post-sawtooth zone indices real*8, intent(out) :: zto(inzmax) ! post-sawtooth material fractions integer, intent(out) :: ierr ! error code (0=OK)