\input texinfo  @c -*-texinfo-*-
@setfilename drat.info
@settitle DRAT Manual
@dircategory Miscellaneous
@direntry
* Drat: (drat).       Yorick-based 2D ray transport post-processor.
@end direntry
@setchapternewpage odd
@comment %**end of header (This is for running Texinfo on a region.)

@titlepage
@sp 6
@center @titlefont{DRAT Manual}
@sp 1
@center A portable ray transport post-processor
@sp 4
@center December 1990
@sp 5
@center David H. Munro
@sp 1
@center Lawrence Livermore National Laboratory
@page
@vskip 0pt plus 1filll
No fine print --- what you see is what you get.
@end titlepage
@page

@node Top, Introduction, (DIR), (DIR)

@menu
* Introduction:: Availability and usage of DRAT.
* Coordinates::  Tracking a 3-D ray through a 2-D cylindrical mesh.
* Transport::    Solving the transport equation for a given ray path.
* Examples::     A few typical cases.
@end menu

@node Introduction, Coordinates, Top, Top
@chapter Availability and Usage of DRAT

@dfn{DRAT} is a post-processor for LASNEX which integrates the radiative
transport equation along arbitrary rays, according to the opacities,
emissivities and geometry saved in Q-dump files.  DRAT is portable; it
is written in ANSI standard C and uses PDB files for both input and
output.  We have successfully compiled and run many such programs on Sun
and DEC workstations, MacIntoshes, and Crays.  A Sun version of DRAT is
available on the X-division network in @file{/home/risky/munro/bin/drat},
and a Cray version is in my take directory @file{.634366:drat:drat}.
DRAT is intended to replace the DIRT and TDG codes.  It should give
identical results on identical problems, modulo roundoff differences
(which can range up to several 100 ppm due to the half precision opacity
data stored in the Q-dump files).

The command line for DRAT is:
@example
drat @var{rayfile} @var{DumpFileQ00} [-o @var{OutputFile}]
@end example

You must run TOPDB on the Q-dump files produced by LASNEX in order to
convert them to the PDB format @var{DumpFileQ00} (the first file of a
family).  TOPDB is @file{/home/risky/munro/bin/topdb} on the Suns or
@file{.634366:pdb:topdb} under XP on the Crays.  (Start TOPDB with no
arguments for an explanation of its use.)  The Q-dump file(s) contain
the following arrays:

@table @asis
@item integer kmax, lmax, ngroup
The dimensions of the mesh and the number of groups.

@item integer isymz, khold, lhold, kxlm
If the problem had reflection symmetry about z=0, isymz=1,
else isymz=0.  If there are hold lines in the problem, they
are at khold and lhold.  The kxlm dimension is kmax*lmax,
except in very old Q-dump files, where it may be kmax*lmax+1.

@item double gb(ngroup+1), gav(ngroup)
The photon group boundary energies are gb (keV);
gav (keV) are the group center energies.

@item double time; integer ncyc
The problem time (sh) and cycle number.  These and the following
variables are repeated once per dump.

@item double zt(kmax,lmax), rt(kmax,lmax)
@itemx integer ireg(kmax,lmax)
The mesh coordinate arrays (cm) and region numbers.

@item double akap(kmax,lmax,ngroup), ekap(kmax,lmax,ngroup)
The opacity array is akap (1/cm), the source function is
ekap (jrk/sh/cm2/ster/keV).
@end table

Additionally, you must use DAP to write a PDB file, the @var{rayfile},
which contains a description of the rays on which you want the transport
equation integrated.  More than one @var{rayfile} may appear on the
command line.  One @var{rayfile} file must contain, at minimum,
variables with the following names:

@table @asis
@item double x(nrays), y(nrays), z(nrays), theta(nrays)
The rays pass through the points (x,y,z) (cm), making angles
theta (radians) from the +z-axis (@pxref{Coordinates}).
These variables may actually have any number of dimensions,
or be scalars, so long as all four are the same.
@end table

Optionally, @var{rayfile} may include variables with any of the following
names:

@table @asis
@item integer zsym, sphere, flate
If present, zsym overrides the isymz flag in the Q-dump files.
The sphere flag should be present and set to 1 if the problem
is a 1-D sphere (there is no sure way to determine this from
the Q-dump files).  If flate=0, the source function is a
piecewise linear function along the ray path (the default
behavior); if flate=1, the source function is constant within
each zone along the ray path.

@item double tlimits(2)
The transport calculation will not begin until a dump with
time>=tlimits(0) is reached, and will stop before processing
a dump with time>tlimits(1).

@item double slimits(2,nrays)
The integration limits for each ray are from minus to plus
infinity by default, or if slimits(0,i)>=slimits(1,i), else
the transport integration begins at slimits(0,i) and ends
at slimits(1,i) (cm) (@pxref{Coordinates}).  If present,
slimits must be in the same file as x, y, z, and theta.

@item blist
A list of 0-origin indices into the gb array used to restrict
DRAT to a subset of the photon groups present in the Q-dump
files.  The blist should be a monotone increasing list of
numbers >=0 and <=ngroup, which increments by one for each
group to be included in the transport calculation, or by
more than one to skip over some groups.

@item double akapm(kmax,lmax,ngrp), ekapm(kmax,lmax,ngrp)
Zone and photon energy dependent opacity (akapm) and source
function (ekapm) multipliers.  The ngrp index is either
ngroup, if no blist is present, or the number of groups
specified by the blist.

@item integer nomilne(2,2,nsegs)
Prevent DRAT from applying a Milne boundary condition on
the specified logical line segments (the nsegs dimension
may be omitted if there is only one such segment).  Each
of the nsegs segments is a 2-by-2 array [[k1,l1],[k2,l2]],
where k1, l1, k2, and l2 are intended to be the arguments
on a NORAD S-card in the LASNEX generator deck.  The
location of the NORAD sources is not in the Q-dump files,
though it should be.  DRAT never applies a Milne condition
if flate=1, so nomilne is ignored in this case.
@end table

DRAT writes its results to a PDB time history file family.  You must use
DAP to plot or print the results.  The input x, y, z, theta, slimits,
gb, gav, and blist variables are copied as static variables to the
output file.  If blist is present, an additional index array glist,
which is the 0-origin indices of the active groups in the gav array, is
also written as a static variable.  The dynamic variables written to
each record of the output file are:

@table @asis
@item double time; integer ncyc
The problem time (sh) and cycle number.

@item double trfrac(ngrp,nrays)
The multigroup transmission fraction for each active group
(ngrp=ngroup if blist is not present), for each ray.

@item double bright(ngrp,nrays)
The self emission (jrk/sh/cm2/ster/keV) for each active group
(ngrp=ngroup if blist is not present), for each ray.
@end table

DRAT is stripped to the essentials.  Frequently you need to add filter
and instrument response, in photon energy, time, or space, or to include
the effects of a backlighter.  Given the full frequency, ray, and time
dependent trfrac and bright arrays in the DRAT output file, performing
these additional calculations with DAP should be straightforward.  Since
DRAT works entirely through PDB files, you can switch between the Crays
and workstations at any point you wish.

@node Coordinates, Transport, Introduction, Top
@chapter Coordinate Systems for Mesh and Rays

DRAT represents a break from the coordinate system used in the older
DIRT and TDG codes.  The new ray coordinate system is easier to
describe, but the general procedure for conversion from DIRT or TDG
coordinates to DRAT coordinates is tricky.

LASNEX works with a quadrilateral mesh in a cylindrical coordinate
system, @i{(z,r)}.  This coordinate system is conventionally plotted
with the @i{+z}-axis increasing to the right, and the @i{r}-axis
increasing upward.  (Negative values of @i{r} are meaningless and do not
occur in LASNEX.)  Each quadrilateral zone actually represents a torus
of material; each edge of the mesh is a segment of a cone.  The 3-D
nature of the mesh is implicit in LASNEX itself, but becomes explicit
for DRAT, which must track an arbitrary straight line through all of
these conical toruses.  To handle this, DRAT enhances the LASNEX by
adding the @i{x} and @i{y} coordinates such that:
@example
x^2 + y^2 = r^2
@end example
@noindent
and @i{(x,y,z)} is a right-handed orthogonal coordinate system.

A ray in DRAT travels in a plane of constant @i{y}, at an angle
@i{theta} from the @i{+z}-axis, so that if @i{theta = +pi/2}, the ray
travels in the @i{+x}-direction.  You specify a ray by giving
@i{(x,y,z)} for any point on the ray, and @i{theta}.  DRAT denotes the
coordinate along the ray by @i{s}, so that a photon starts at @i{s =
-infintity} and flies to @i{s = +infinity}.  The origin @i{s = 0} is
where the plane perpendicular to the ray contains the origin of the
@i{(x,y,z)} coordinate system.  (Formally, @i{s = z*cos(theta) +
x*sin(theta)}.)

The most convenient way to think about the DRAT coordinate system is to
imagine that the @i{+x}-direction coincides with the @i{r}-axis of a
standard LASNEX mesh plot.  The @i{+y}-axis is out of the paper toward
your face, and DRAT rays travel in planes of constant @i{y}, in a
direction @i{theta} radians counterclockwise from the @i{+z}-direction.

The @i{s} coordinate along the ray comes into play if and only if you
specify an slimits variable in your @var{rayfile}.  If slimits is
present, DRAT starts integrating the transport equation at @i{s} =
slimits(0), and stops at @i{s} = slimits(1) (you must supply one
slimits pair for each ray).  Unlike tlimits (the problem time limits),
which have the granularity of the dump times, DRAT applies slimits
precisely, interpolating to the interior of zones if necessary.

@node Transport, Examples, Coordinates, Top
@chapter Integrating the Transport Equation

The Q-dump files contain the kmax-by-lmax mesh arrays @kbd{zt} and
@kbd{rt}, as well as the kmax-by-lmax-by-ngroup opacity array
@kbd{akap} and emissivity array @kbd{ekap}.  The units of @kbd{zt}
and @kbd{rt} are cm, @kbd{akap} is in 1/cm, and @kbd{ekap} is in
jrk/sh/cm2/ster/keV (that is, specific intensity units, the same as
sigma times temperature to the fourth power).  The differential
equation for the specific intensity @kbd{Inu} along a ray is:

@example
 dInu
------- = -Inu + ekap
akap*ds
@end example

This equation can be more conveniently expressed for numerical purposes
as a recursion relation; each chunk of material transmits a certain
fraction @kbd{trfrac} of the radiation incident on it along the ray,
and contributes some self-emission @kbd{bright} as the ray emerges. If
chunk 2 follows chunk 1, then:

@example
bright(1 then 2) = bright(2) + trfrac(1)*bright(1)
trfrac(1 then 2) = trfrac(2)*trfrac(1)
@end example

In DRAT, the chunks are the segments of the ray within each zone along
the ray path.  For each zone, the transmitted fraction is given by
@example
trfrac = exp( -akap(zone) * ds(zone) )
@end example
@noindent
where @kbd{ds} is the length of the segment of the ray lying within the
zone.  (Note that a ray is a straight line in 3-D, while the edges of a
zone are sections of cones.)  That is, DRAT assumes that the opacity is
a step function along the ray, assuming its zone centered value at all
points of the ray lying within a zone.  For the self-emission from a
zone, you have a choice of two models, which you select using the
@kbd{flate} switch (@pxref{Introduction}).

If @kbd{flate=1}, DRAT assumes that the source function @kbd{ekap} is
also a step function along the ray, with its zone centered value at all
points of the ray within each zone.  The self-emission from a zone is
then simply:
@example
bright = (1 - exp(-akap(zone)*ds(zone)) * ekap(zone)
@end example

By default (or @kbd{flate=0}), DRAT point centers the @kbd{ekap} array
computed by LASNEX, then linearly interpolates @kbd{ekap} along each
edge to the point where the ray cuts the edge.  In this mode, DRAT
assumes that the source function varies linearly along the ray path from
the value as it cuts the edge entering a zone to the value as it cuts
the edge leaving the zone.  Solving the transport equation for
@kbd{ekap} varying linearly from @kbd{ek1} to @kbd{ek2} over the length
@kbd{ds} with @kbd{akap} constant yields:
@example
dtau = akap(zone) * ds(zone)
bright = (1-exp(-dtau))*ek1 + (1 - (1-exp(-dtau))/dtau)*(ek2-ek1)
@end example

DRAT uses an extremely computation intensive point centering algorithm
copied from the TDG and DIRT codes.  The value of the point centered
@kbd{ekap} at a point is the weighted average of the (up to four) zones
surrounding the point.  The weighting factor is:
@example
((1-exp(-akap*ddk))^2 + (1-exp(-akap*ddl))^2) / (akap*area)
@end example
@noindent
where @kbd{area} is the zonal area and @kbd{ddk} and @kbd{ddl} are the
median lengths of the zone in the k and l directions.  This function
weights zones with optical depths near unity most heavily.  For points
along a non-symmetry problem boundary which is not a hold line, DRAT
applies an additional Milne boundary condition, again following TDG and
DIRT.  This condition is certainly incorrect if the boundary had a NORAD
source in the LASNEX problem, since the Milne condition only makes sense
for a vacuum boundary.  Unlike TDG or DIRT, DRAT allows you to correct
this mistake by specifying the NORAD boundary using the @kbd{nomilne}
variable (@pxref{Introduction}).

The Milne condition reduces the source function at the boundary points
if the zone is optically thick in the direction (logically) perpendicular
to the boundary.  If the optical depth perpendicular to the boundary is
@kbd{dtau}, then the source function @kbd{ekap} is reduced by the factor
@kbd{1+1.5*dtau}.  If @kbd{dtau>1}, DRAT checks the gradient of @kbd{ekap}
(again logically perpendicular to the boundary), and may reduce the
@kbd{dtau} (but never to less than 1) if the gradient is small (i.e.- if
LASNEX does not seem to be losing much energy out the boundary).

The point centering scheme and especially the Milne condition are of
rather dubious pedigree.  Besides being expensive and probably wrong,
the default @kbd{flate=0} model also gives results which are
discontinuous with respect to small changes in ray position.  Curiously
enough, the simple @kbd{flate=1} model with its step function source
gives results which vary continuously with ray position.  The resolution
of this seeming paradox is that when a ray just begins to touch some
zone, a new interpolation point is added to the piecewise linear source
function, which discontinuously changes it.  This happens even with the
simplest ray geometries, such as 1-D spheres.  If no one can think of a
good reason other than historical inertia, I may well decide to make
@kbd{flate=1} the default.  Until then, I encourage you to specify
@kbd{flate=1}.

@node Examples, , Transport, Top
@chapter Examples

My own favorite problems are Rayleigh-Taylor simulations.  I run LASNEX
with an 11-by-38 mesh (418 points, 370 zones), which is a torus at very
large @i{r}, representing half a wavelength of a perturbation in the
@i{r} direction, which is roughly the k direction.  The @i{z} direction
is the direction of motion of the ablating foil.  The k= 1 line is a
hold line, and I employ a NORAD source (and a VEL source) at k= 11.  I
typically want to integrate the transport equation for 10 rays parallel
to the @i{z}-axis, with one ray centered in each column of zones.  If
@i{r}= 1.0025 cm at k= 1 and @i{r}= 1.0000 cm at k= 11, then I can use
DAP to set up several @var{rayfiles} as follows:

@example
miggle[5] printenv BASIS_CODE_DIR
/usr/local/vbasis/lib/basis
miggle[6] xdap
DAP> double x(10), y(10)=0, z(10)=0, theta(10)=0
DAP> x= 1.000125 + 0.00025*(iota(10)-1)
DAP> integer nomilne(2,2)= [[11,1],[11,37]]
DAP> osetfile rays.pdb, pdb
DAP> store x,y,z,theta,nomilne
DAP> integer flate=1
DAP> osetfile flate.pdb, pdb
DAP> store flate
DAP> end

   CPU (sec)    SYS (sec)
       1.140       0.250
miggle[7] 
@end example

The file @file{rays.pdb} contains the ray coordinates and specifies the
NORAD boundary condition on k= 11 (the hold line at k= 1 is specified in
the Q-dump files).  The file @file{flate.pdb} simply contains a flate= 1
switch.  Note the environment variable @kbd{BASIS_CODE_DIR}, which must
be set in order for DAP to work.  Also, DAP requires you to declare
floating point numbers as ``double'' on the workstations, but ``real''
on the Crays.  If the Q-dump files begin with @file{rth79Q01}, then I
run DRAT with:

@example
miggle[11] drat rays.pdb rth79Q01 -o rth79Q01.l
38.6u 1.6s 0:40 98% 0+976k 1+25io 0pf+0w
miggle[12] drat rays.pdb flate.pdb rth79Q01 -o rth79Q01.f
8.6u 1.5s 0:10 95% 0+936k 6+27io 0pf+0w
miggle[13] 
@end example

The results of the first run will be in @file{rth79Q01.l}.  The results
of the second run differ by using the flat emissivity model, and are in
@file{rth79Q01.f}.  Notice the large difference in the amount of CPU
time required for the two models (39 seconds as compared to 9 seconds,
on my Sun-4/330).  To see the results, use DAP again:

@example
miggle[13] xdap
DAP> ps rel; win on
DAP> load rth79Q01.l; dbrbase="rec"
DAP> listf
nomilne  gb       theta    gav      slimits
Base record
rec                                             nrecords =    17
        time      ncyc     trfrac   bright
DAP> plot bright
DAP> time
time     =    1.40491D-01
DAP> setrec 11
DAP> plot bright
DAP> time
time     =    2.40155D-01
DAP> load rth79Q01.f; dbrbase="rec"
DAP> setrec 11
DAP> plot bright
DAP> bright(15,)
bright(15,)        shape: (1,10) 
  1:         2.06608D-05   2.09432D-05   2.06591D-05   2.04980D-05
  5:         2.03434D-05   2.04377D-05   2.00107D-05   2.09656D-05
  9:         2.07852D-05   2.15009D-05 
DAP> prtfile
CURRENT INPUT FILE:  rth79Q01.f 
CURRENT RECORD NUMBER:     11
CURRENT RECORD BASE:  rec                     
DAP> setfile rth79Q01.l
DAP> setrec 11
DAP> bright(15,)
bright(15,)        shape: (1,10) 
  1:         2.06927D-05   2.08857D-05   2.07160D-05   2.05281D-05
  5:         2.03687D-05   2.04197D-05   1.99121D-05   2.08169D-05
  9:         2.09583D-05   2.14623D-05 
DAP> setfile rth79Q01.f
DAP> time
time     =    1.40491D-01
DAP> end

   CPU (sec)    SYS (sec)
       1.900       0.630
miggle[14] 
@end example

Note the use of the DAP functions and variables @kbd{load}, @kbd{dbrbase},
@kbd{setrec}, @kbd{setfile}, @kbd{listf}, and @kbd{prtfile}.  In this case,
the problem is fairly well resolved, and the difference between the flat
source function model and the linear source function model is quite
small.

@ignore
@node Index,, Cap Index, Top
@unnumbered Concept Index

@printindex cp
@end ignore

@contents
@bye
