A brief introduction to the RPLOT calculator -- D. McCune 28 Sept 1999. Abstract: ========= The RPLOT calculator is a tool for arithmetic manipulation of results of TRANSP simulations. The calculator permits * formation of arithmetic expressions involving data contained in the TRANSP output database: scalar and profile functions of time. * creation of user defined profile or scalar functions which can be referred to in subsequent expressions, using assign statements or commands. * definition of multigraphs and manipulation of multigraph contents, using commands. The calculator interprets user input, read in the form of a character string (max length 512 characters). It keeps track of the results of the most recent calculation in an "accumulator". However, by using user defined functions, calculations involving many steps can easily be constructed. Basic Syntax Options: ===================== RPLOT calculator inputs take one of the four following forms: (1) a bare expression: for example: 1.5*ne*te -- the arithmetic expression is evaluated point by point; the result is stored in the calculator "accumulator". Since the operands in this expression are "zone centered" profile functions of time and radial flux coordinate, the result of the expression evaluation is also a zone centered profile function of time and radial flux coordinate (as one might expect). If the calculator is being called from a program, the results are returned in the data buffer array provided by the caller. (2) an assign statement using an expression: for example: PTHE = 1.5*ne*te --or, with labeling information: PTHE,"electron thermal energy density","ptcl-ev" = 1.5*ne*te -- the expression is evaluated and the result is stored in the user defined function "pthe". The name must not conflict with any permanent function name or multigraph name. Reuse of user defined function names will cause replacement of data, with a warning. -- optionally, labeling information can be provided. The first label is a description label, the second a physical units label; always enclose in quotes to prevent mistakes by the calculator input parser. The description label can be up to 32 characters long; the units label up to 16 characters long. -- all function names are folded to uppercase. (3) a command: for example: %SAVE(pthe,"electron thermal energy density","ptcl-ev",1.5*ne*te) --or, %MAXPRO(1.5*ne*te) -- the first command saves the results of the expression "1.5*ne*te" in the user defined profile function "PTHE" with the specified description and units label; i.e. equivalent to the 2nd assign statement example in (2) above. -- the second command returns a scalar function of time which is the max of "1.5*ne*te" profile at each time. (4) an assign statement using a command: for example: PTHE_MAX = %MAXPRO(1.5*ne*te) -- this stores %MAXPRO(1.5*ne*te) as a user function with default labeling. Note that only commands with numeric output can be used in this manner. Some commands do not produce numeric output. When writing scripts involving calculator expressions, the character "#" may be used as a comment character, anywhere in the line not inside quotes. In arithmetic expressions, the symbol "$" refers to the current contents of the accumulator. Arithmetic Expressions: ======================= The heart of the RPLOT calculator functionality is the arithmetic expression evaluator. This is built on a parser which recognizes operators and operands and certain punctuation such as parentheses and commas. Operators and operands are discussed separately, operators first. RPLOT calculator operators have one of four syntaxes: MONADIC, POSITIONAL: +, - example: -ne gives the additive inverse of the profile function NE (electron density). DYADIC, POSITIONAL: algebraic operators: partial list: operator description Precedence ** exponentiation 5 * multiplication 3 / division 3 + addition 2 - subtraction 2 example of use: (1.5*(ne*te + ni*ti))**0.5 parentheses are used to modify the precedence order of evaluation. Equal precedence operators are parsed left to right, except chained "**" operators which are parsed right to left. MONADIC, FUNCTIONAL: partial list: SQRT, ABS, VOLINT, GRAD example: VOLINT(ne) (volume integral of electron density, yielding total number of enclosed electrons as a function of flux surface, and time). Some monadic functional operators are applicable to any profile, others are applicable only to an appropriate type of profile, e.g. a function of flux zone and time but not a function of poloidal angle and time. Many operators make implicit reference to metric information stored in the TRANSP database. For example, GRAD is evaluated using information equivalent to <1/dr> which was numerically flux surface averaged and then stored as a named output in the TRANSP database. The complete list of available operators is available by "help" calls from programs which make use of the RPLOT calculator. DYADIC, FUNCTIONAL: partial list: MIN, MAX example: max(ne*te,ni*ti) Evaluate a point by point maximum. Both arguments must be of the same type (i.e. both scalar functions, or, both profile functions vs. time and a particular x axis coordinate). A preparse converts expressions like max(a,b,c) to max(a,max(b,c)), so that, in the case of MIN and MAX only, a multi-adic syntax is allowed. ----- In RPLOT calculator expressions, OPERANDS can be any of the following: ->Constants: a floating point number, or "PI", or "EE" (=exp(1.0)). ->Scalar functions of time: a data function of time only, e.g. PCUR in most TRANSP runs. ->Profile functions of time: a data function of time and some other coordinate, such as a radial flux coordinate (zone or bdy centered). ->"$" refers to the "RPLOT ACCUMULATOR" which contains the result of the most recent expression evaluation. It can be referred to in the next expression. However, evaluation of any expression overwrites "$". ->Any valid (sub)expression can be an operand. There are some restrictions on operands. Generally, dyadic operators require profile operands to be of the same type. If one operand is a profile and the other a scalar or constant, the scalar or constant is "promoted" to be a profile of the same type. The output of an operator is always of the same type as the input operands. In the case of profiles, "same type" means being defined against the same non-temporal coordinate. For certain types of profiles, these rules are relaxed and automatic type conversion is provided. The following automatic type conversions are supported: "zone centered" <--> "bdy centered" "zone centered" --> "midplane major radius" "bdy centered" --> "midplane major radius" Where: "zone centered" refers to a quantity which is associated with the *center* of numerical zones in the TRANSP simulation, "bdy centered" refers to a quantity which is associated with the *outer boundary* of numerical zones in the TRANSP simulation, and "midplane major radius" refers to a quantity which is a function of midplane radius - flux surface intersection locations (or a generalized definition thereof in the case of updown asymmetry). The following diagram may help. The intersection of the TRANSP model plasma with the midplane may be thought of as looking like: | * | * | * | * x * | * | * | * | Rmajor-----> where "|" denotes zone boundaries, "*" denotes zone centers, and "x" denotes the magnetic axis. Except at the axis, each flux surface intersects the midplane twice, on the inboard and outboard side of the plasma model. In this picture, a "zone centered" function is defined at the points {"*"}, and is presumed to be flux zone constant, i.e. the inboard and outboard values or the same, and the number of points at each time point is four. (Real TRANSP runs generally have 20 or 50 zone centers). A "boundary centered" function also contains four points, covers the points {"|"} but not the axis point "x". A function of major radius (in this example picture) would contain 9 points and would cover the boundaries on both sides plus the magnetic axis. (Real TRANSP runs have midplane major radius grids of 41 or 101 points). The automatic transformation between these coordinate grids involves interpolation, and usually also a half zone width extrapolation at at (at least one of) the edges; because of this, the calculator will write an informational message about such transformations. The back transformation from major radius to flux surface or flux zone is not automatic, because, generally, functions of major radius are not flux surface constant, and so, the user needs to specify a "symmetrization" method; see %XMAP under the section on commands. Expression Evaluation and Error Trapping: ========================================= The Expression Evaluator works by carrying out a stack based "reverse polish" parse of the expression. This leads to an easily evaluated stack based method for carrying out the computation in the order as specified either by precedence rules or explicit use of parentheses in the input expression. Expressions are evaluated in a loop over time points, and each operator returns a result of the same type as its operands. This means that operators which need to refer across time points (such as time integration or time differentiation), and operators which need to have a different type of output than its inputs, are not supported. Such operators are implemented as calculator "commands" instead (see below). During evaluation, operations (such as SQRT, exponentiation, and divide) that can lead to arithmetic errors are checked. Thus, arithmetic errors are trapped with a warning message to the user. When arithmetic errors occur, only portions of the resulting calculation are valid. Generally, users should fix such situations e.g. with MIN or MAX functions. For example: BT = FBTX*BZXR/RMAJM # |Bt| vs Rmajm BTOT = FBX*BZXR/RMAJM # |B| vs Rmajm BP = SQRT(max(0.0,(BTOT**2-BT**2))) # |Bp| vs Rmajm The MAX would prevent an error due to SQRT(negative) in case of round of error during evaluation when B=~BT. Note: a better way to calculate BP on the midplane would be to use TRANSP's FBPBT function: BP = FBPBT*FBTX*BZXR/RMAJM # here BP is signed, though. In these expressions, RMAJM is the major radius grid, BZXR is the external toroidal field in its scalar form R*BT, and FBX, FBTX, and FBPBT are, respectively, |B|/|BT(ext.)|, |BT|/|BT(ext.)|, and |BP|/|BT|, all functions of midplane major radius defined by TRANSP. One could also try TMP=%RMJMAP(PLFLX) # create poloidal flux, webers, vs. Rmajor, cm BP = (1.0E4/RMAJM)*D(PLFLX)/D(RMAJM) # note units conversion factor 1.0e4 yielding the poloidal field, tesla, vs. major radius. (The generic finite difference operator, D(...), is new). Calculator Commands: ==================== Calculator "commands" are used to implement * arithmetic operations not compatible with the arithmetic expression evaluation scheme * "meta" operations such as creation/deletion of user defined functions and multigraphs, or access to data from a different TRANSP run. Although some commands accept expressions as arguments, even commands which have arithmetic output cannot be used directly in expressions. Thus: TMP1=%MAXPRO(ne*te) TMP2=%MAXPRO(ni*ti) FUNC=MAX(TMP1,TMP2) but not FUNC=MAX(%MAXPRO(ne*te),%MAXPRO(ni*ti)) # ?? ** invalid syntax ** "Meta" commands, which have no arithmetic output, cannot appear on the right hand side of an assign statement. Thus: %DELETE(*) # delete all user defined functions but not TMP=%DELETE(*) # ?? ** invalid, %DELETE has no output. All commands have arguments, which can be specified either positionally or by keyword. Most arguments have default values. Generally, defaulted arguments can be omitted. In the case of positional syntax, the correct number of commas must be supplied up to the last explicitly given argument. For example, the %SAVE command takes four arguments with the given ordering: FCN_ID,LABEL,UNITS,EXPR. Thus, %SAVE(FOO) # save accumulator as user variable FOO, default labeling. %SAVE(FOO,,,$) # same result, accumulator explicitly named. %SAVE(FCN_ID=FOO) # same result, keyword syntax %SAVE(FCN_ID=FOO,EXPR=$) # same result, $ explicitly named. But %SAVE(FOO,$) # Save accumulator with description label "$" means something different and probably not intended. The following is a synopsis of available commands. "EXPR" always means arithmetic expression; any arbitrarily complicated RPLOT calculator arithmetic expression may be used. If a command has EXPR as an argument, it has as numeric output either EXPR itself or a transformation thereof, and the accumulator is updated. The RFETCH command also updates the accumulator with the data it reads from a separate TRANSP run results archive. SUMMARY of COMMANDS: -> creation / deletion / relabeling of user functions %SAVE(FCN_ID,LABEL,UNITS,EXPR) ... create a user defined function %DELETE(FCN_ID) ... delete a user defined function %LABEL(ITEM_ID,LABEL,UNITS) ... relabel existing function or multigraph -> creation / deletion / modification / relabeling of multigraphs %MG_CREATE(PKG_ID,LABEL,FCN1,FCN2,FCN3,FCN4,FCN5,FCN6) ... create a multigraphs with label and ... up to 6 functions, FCN1 must be given. ... all FCNs must be of the same type, ... all must have the same physical units. ... defaulted FCN arguments are ignored. %MG_ADDFUN(PKG_ID,FCN1,FCN2,FCN3,FCN4,FCN5,FCN6,FCN7) ... add FCNs to an existing multigraph. %MG_DELFUN(PKG_ID,FCN1,FCN2,FCN3,FCN4,FCN5,FCN6,FCN7) ... remove FCNs from an existing multigraph. %MG_DELETE(PKG_ID) ... delete a user defined multigraph %LABEL(ITEM_ID,LABEL,UNITS) ... modify label of a multigraph. -> arithmetic "reductions" (profile -> scalar function of time): %TIME_TRACE(INDEX,EXPR) ... time variation at INDEX'th x point. %XINTERP(X,EXCEPTION,EXPR) ... interpolation to (time dependent) ... X location; value=EXCEPTION when ... X is not in range (Xmin,Xmax). %XLOCATE(TEST_VALUE,NOT_UNIQUE,NOT_FOUND,EXPR) ... return as function of time the X ... location of TEST_VALUE in the ... profile EXPR; else return an ... exception number NOT_UNIQUE or ... NOT_FOUND. TEST_VALUE can be ... set to MIN or MAX or a numeric ... value (example: XQ32=%XLOCATE(1.5,,,Q) TEQ32,"Te at q=3/2","EV" = %XINTERP(XQ32,,TE) returns electron temperature at the q=1.5 surface vs. time. A default value of -99.0 is returned at those times when the q=1.5 surface is either non-existant or non-unique). %MINPRO(EXPR) ... return MIN(profile) as function of time %MAXPRO(EXPR) ... return MAX(profile) as function of time -> arithmetic operations involving time: %TIME_DERIV(EXPR) ... numerically evaluate d/dt(EXPR) %TIME_INT(T0,EXPR) ... numerically integrate from t=T0 ... result is function of time t' = the ... upper limit of the integration. %SMOOTH(DELTA_T,DELTA_X,EPS_T,EPS_X,EXPR) ... epsilon/delta triangular smooth of EXPR ... DELTA_X and EPS_X ignored when EXPR ... is a function of time only. %TIME_AVG(DELTA_T,EXPR) ... +/- DELTA_T point by point time average ... of EXPR. -> mapping (type change) operations: %XMAP(TARGET,SIDE,EXPR) ... map function of major radius to flux ... surface; TARGET=CTR or BDY; ... SIDE=IN,OUT or INOUT for averaging. %RMJMAP(EXPR) ... map flux surface/zone function to ... major radius -> fetch data from a different TRANSP run: %RFETCH(PATH,RUN_ID,FCN_ID,LOCAL_ID) ... open TRANSP run (RUN_ID) at location ... (PATH) and read that run's function ... (FCN_ID). Interpolate to current ... run's timebase and store as user ... defined function LOCAL_ID. (example #1a: %RFETCH("/u/transp/result/TFTR.88",37065Z15,PCUR,PCUR2) means go to the indicated TRANSP run in the indicated directory, read out the function PCUR, interpolate this to the current run's timebase, and call the result PCUR2; quotes suppress uppercase conversion). (example #1b: %RFETCH(".",37065Z15,PCUR,PCUR2) means fetch from the same source (possibly MDS+) where the session's main run data was read). (example #1c: %RFETCH(" ",37065Z15,PCUR,PCUR2) means read from the current working directory). (example #1d: %RFETCH("./work",37065Z15,PCUR,PCUR2) means read from the work subdirectory of the current working directory, on unix systems). (example #1e: %RFETCH(".TFTR.88",37065Z15,PCUR,PCUR2) means read from $RESULTDIR/TFTR.88 (unix) or RUNDATA:[TRANSP.TFTR.88] (vms) in the current file system-- following the old RPLOT shorthand). (example #2: %RFETCH("MDS+:TRANSPGRID.PPPL.GOV:TRANSP:JET.00",50632C02,TE,TE2) means go to the indicated TRANSP run on the indicated MDS+ server, read out the function TE, interpolate this to the current run's timebase, and call the result TE2. The syntax is for PPPL-style TRANSP tree identification, with treename=TRANSP always, and with tok.yy and runid, instead of shot number). (example #3: %RFETCH("MDS+:CMODA.PSFC.MIT.EDU:TRANSP03",12345678,NE,NE2) means go to the indicated TRANSP run on the indicated MDS+ server, read out the function NE, interpolate this to the current run's timebase, and call the result NE2. The syntax is for CMOD-style TRANSP tree identification, using shot number and different treenames for different runs on the same shot). Note that the first field after "MDS+:" is a server name; it is allowed to also specify a port, as in "EUROPA.PPPL.GOV:8501". In this case, the port number counts as part of the server name specification. The general RFETCH MDS+ argument syntax, then, is MDS+:::, , , or MDS+::, , , ---------------------> new September 2001 <------------------------- -> fetch run from a GS2 tabulated file (created by R. Budny method): %GS2FETCH(PATH,XID,PREFIX) PATH gives a file or MDS+ path to GS2 tabulated results data XID choose selects flux coordinate "X" or major radius "RMAJM" PREFIX (max 5 characters) determines the names of three profile functions created as a result of successful completion of the command: _AKY -- normalized k(theta) _OMEGA -- normalized frequency _GAMMA -- normalized growth rate the precise meanings of these items should be obtained via consultation with an expert. Concluding Remarks. =================== This document describes an upgraded version of the RPLOT calculator which has been used for years by experimental physicists as a tool for examining results of TRANSP simulations of tokamak experiments. The new calculator is accessible via the traditional (command line oriented) "rplot" graphics program, but a callable interface is also available. The description of this interface is in the documentation of the "trread" module at the NTCC modules library website: -- go to http://w3.pppl.gov/ntcc -- click on "Modules Library" -- click on "Data Analysis and Visualization" -- click on "trread" -- click on the link to the trread module documentation. The hope is that this callable interface, from which all dependence on legacy command line interface or PPPL native tektronix graphics display routines has been removed, will be used as a foundation for a modern GUI based TRANSP results visualization tool. Almost all of the functionality found in the traditional command line "rplot" program can be accessed using this interface (exceptions are: SNAP and UFILES i/o features). The calculator supports a method for creating a scripted series of arithmetic calculations. It is not a "full fledged programming language" for a number of reasons, most notably being that it lacks constructs for looping and conditional branching. The author feels that the effort required to develope and maintain a full fledged programming language would not be warranted. Questions / Comments should be emailed to the author: Doug McCune dmccune@pppl.gov