Google

i/test1

yorick banner

Home

Manual

Packages

Global Index

Keywords

Quick Reference


/*
   TEST1.I
   Poorly vectorizing physics problem for timing Yorick.

   $Id: test1.i,v 1.1 1993/08/27 18:50:06 munro Exp $
 */
/*    Copyright (c) 1994.  The Regents of the University of California.
                    All rights reserved.  */

/* WARNING -- this code assumes 1-origin indexing is the default */


func test1 (npass)
/* DOCUMENT test1
         or test1, npass
     Track a mock "ablation front" as it propagates through a mesh.
     If NPASS is given, the calculation is repeated that many
     times.  The zoning, densities, temperatures, pressures, and
     velocities are all computed arbitrarily, but the number of zones
     and groups are taken to be representative of a typical 1-D
     ablation calculation.  */
{

  if (is_void(npass) || npass<=0) npass= 1;

  now= split= array(0.0, 3);
  timer, now;
  n= npass;
  while (n--) {
    prevScale= prevAbl= 37;

    for (i=1 ; i<=ntimes ; i++) accum, i;
  }
  timer, now, split;

  timer_print, "Time per pass", split/npass, "Total time", split;

  if (prevAbl!=15 || prevScale!=8 ||
      nallof(approx_eq(rorab,
		       [2.673480885e-03,2.566922597e-03,2.422313787e-03,
			2.211817497e-03,1.870390096e-03])) ||
      nallof(approx_eq(times,1.5)) ||

      nallof(approx_eq(rordx,3.109107238e-03)) ||

      nallof(approx_eq(rorsx,2.946068988e-03)) ||

      nallof(approx_eq(rorhx,2.916081687e-03)) ||

      nallof(approx_eq(dmax,3.034383698e-01)) ||

      nallof(approx_eq(zs,5.973585920e-03)) ||

      nallof(approx_eq(vel,-1.662706089e-01)) ||

      nallof(approx_eq(acc,-3.148890163e-03)) ||

      nallof(approx_eq(scl,2.061472601e-03)))

    write, "***WARNING*** values returned by accum are not correct";
}

func approx_eq(x, y)
{

  return (abs(x-y)/(abs(x+y)+1.e-35))<1.e-6;
}

#if 0
CPU seconds per pass:
            IDL    Yorick     Basis       DAP
HP730        -      0.41       8.21      8.15
Solbourne  0.57     1.00       20.1      22.4  (1st trial)
Solbourne  0.57     0.90       19.5      20.9  (2nd trial)
Note:  Times on Solbourne fluctuate by ~10%

Run on forum (Solbourne) 8/Dec/92
     forum[4] time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
     IDL. Version 2.3.0 (sunos sparc).
     Copyright 1989-1992, Research Systems, Inc.
     All rights reserved.  Unauthorized reproduction prohibited.
     Site: 1491.
     Licensed for use by: LLNL - X Division

     IDL> .rnew /home/miggle/munro/Yorick/include/test1
     % Compiled module: PEAK.
     % Compiled module: DERIV.
     % Compiled module: ACCUM.
     % Compiled module: APPROX_EQ.
     % Compiled module: TEST1.
     % Compiled module: INTERP.
     % Compiled module: UNCP38.
     % Compiled module: $MAIN$.
     IDL> test1,1
     IDL> exit
     1.1u 0.3s 0:28 5% 0+1264k 4+0io 5pf+0w
     forum[5] !!
     time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
     IDL. Version 2.3.0 (sunos sparc).
     Copyright 1989-1992, Research Systems, Inc.
     All rights reserved.  Unauthorized reproduction prohibited.
     Site: 1491.
     Licensed for use by: LLNL - X Division

     IDL> .rnew /home/miggle/munro/Yorick/include/test1
     % Compiled module: PEAK.
     % Compiled module: DERIV.
     % Compiled module: ACCUM.
     % Compiled module: APPROX_EQ.
     % Compiled module: TEST1.
     % Compiled module: INTERP.
     % Compiled module: UNCP38.
     % Compiled module: $MAIN$.
     IDL> test1,21
     IDL> exit
     12.5u 0.2s 0:30 42% 0+1432k 0+0io 0pf+0w


     forum[6] yorick

      Yorick ready.  For help type 'help'

     > #include "/home/miggle/munro/Yorick/include/test1.i"
     > test1
		    Timing Category     CPU sec  System sec    Wall sec
		      Time per pass       1.010       0.000       1.020
			 Total time       1.010       0.000       1.020
      -----Total Elapsed Times-----       1.670       0.140      20.270
     > test1,21
		    Timing Category     CPU sec  System sec    Wall sec
		      Time per pass       1.001       0.001       1.004
			 Total time      21.020       0.030      21.080
      -----Total Elapsed Times-----      22.710       0.180      51.760
     > test1
		    Timing Category     CPU sec  System sec    Wall sec
		      Time per pass       0.990       0.000       0.990
			 Total time       0.990       0.000       0.990
      -----Total Elapsed Times-----      23.720       0.190     163.440
     > test1
		    Timing Category     CPU sec  System sec    Wall sec
		      Time per pass       0.990       0.010       1.000
			 Total time       0.990       0.010       1.000
      -----Total Elapsed Times-----      24.740       0.200     171.340
     > quit
     24.7u 0.2s 3:35 11% 0+1088k 18+0io 64pf+0w


     forum[7] basis
     Basis    (basis, Version 921125)
     Run at 10:06:44 on 12/08/92 on the forum    machine, suffix 10470x
     Initializing Basis System
     Basis 7.0
     Initializing EZCURVE/NCAR Graphics
     Ezcurve/NCAR 2.0
     Basis> echo=0

     Basis> read /home/miggle/munro/Yorick/include/test1.bas

     End of input file /home/miggle/munro/Yorick/include/test1.bas
     Resuming input from TERMINAL
     Basis> test1(1)

	CPU (sec)   SYS (sec)
	   20.150       0.000
     Basis> end

	CPU (sec)   SYS (sec)
	   23.117       0.467
     23.1u 0.4s 1:12 32% 0+2408k 20+0io 229pf+0w


     forum[8] dap

     DAP> read /home/miggle/munro/Yorick/include/test1.bas
     DAP> test1(1)

	CPU (sec)   SYS (sec)
	   22.383       0.017
     DAP> end
     26.7u 1.1s 1:07 41% 0+2784k 17+0io 209pf+0w


Run on tonto (HP730) 7/Dec/92
     tonto[8] basis
     Basis    (basis, Version 921125)
     Run at 16:06:31 on 12/07/92 on the tonto    machine, suffix 2545x
     Initializing Basis System
     Basis 7.0
     Initializing EZCURVE/NCAR Graphics
     Ezcurve/NCAR 2.0
     Basis> echo=0

     Basis> read /home/miggle/munro/Yorick/include/test1.bas

     End of input file /home/miggle/munro/Yorick/include/test1.bas
     Resuming input from TERMINAL
     Basis> test1(1)

	CPU (sec)   SYS (sec)
	    8.260        .010
     Basis> test1(10)

	CPU (sec)   SYS (sec)
	   82.090        .020
     91.7u 0.1s 4:07 37%


     tonto[9] dap
     DAP> echo=0

     DAP> read /home/miggle/munro/Yorick/include/test1.bas
     DAP> test1(1)

	CPU (sec)   SYS (sec)
	    8.130        .000
     DAP> test1(10)

	CPU (sec)   SYS (sec)
	   81.520        .050
     DAP> end
     91.5u 0.6s 3:02 50%


     tonto[10] yorick

      Yorick ready.  For help type 'help'

     > #include "/home/miggle/munro/Yorick/include/test1.i"
     > test1
		    Timing Category     CPU sec  System sec    Wall sec
		      Time per pass       0.400       0.000       0.410
			 Total time       0.400       0.000       0.410
      -----Total Elapsed Times-----       0.720       0.070      14.810
     > test1,20
		    Timing Category     CPU sec  System sec    Wall sec
		      Time per pass       0.412       0.000       0.599
			 Total time       8.240       0.000      11.980
      -----Total Elapsed Times-----       8.960       0.070      32.350
#endif

/* ------------------- Ablation front tracker ---------------------- */

/* The ablation front location (rorab) is defined here by the condition that
   the matter temperature reaches a particular fraction of the current drive
   temperature.  Five different hypotheses about this fraction (.4 thru .8)
   are calculated simultaneously.  These can be compared with the ablation

   front location (rorhx) defined by the location of peak phot (X-ray
   absorption), or defined by the steepest gradient (rorsx), or to the
   location of maximum density (rordx).  Other useful parameters are also
   recorded for later perusal.  */

func accum (irec)
{

  /* Input arrays are external -- they would ordinarily be read from
     successive time history dumps of a hydro program.  */
  /* time-independent input quantities */
  extern kmax, lmax;   /* number of points for mesh arrays */
  extern fraction;     /* list of fractions for front definition */
  extern ror;          /* column density coordinate */

  /* Since the direction of the ablation is known, the various searches

     begin at the location of the fron on the previous call to accum;
     these are stored in the following two external variables, which

     must be initialized to 37 before the loop calling accum begins.  */
  extern prevScale, prevAbl;

  /* time-dependent input quantities */
  extern time;         /* current problem time */

  extern zt;           /* position (point centered mesh array) */

  extern v;            /* velocity (point centered mesh array) */

  extern s;            /* radiation temperature (zone centered mesh array) */

  extern e;            /* matter temperature (zone centered mesh array) */

  extern d;            /* density (zone centered mesh array) */

  extern p;            /* pressure (zone centered mesh array) */

  extern h;            /* heating (zone centered mesh array) */

  /* output quantities */
  extern times;        /* list of problem times */
  extern rorab;        /* nfraction-by-ntimes ablation depths (rho*r) */

  extern rordx;        /* ntimes rho*r at peak density */
  extern rorsx;        /* ntimes rho*r at steepest density gradient */

  extern rorhx;        /* ntimes rho*r at peak phot (heating rate) */

  extern dmax;         /* ntimes peak densities */
  extern zs;           /* nfraction-by-ntimes ablation front positions */
  extern vel;          /* ntimes unablated slab velocities */
  extern acc;          /* ntimes unablated slab accelerations */
  extern scl;          /* ntimes density scale lengths at front */

  mid = (kmax+1)/2;    /* slab is ablating in l-direction -- calculate
			  front for k-zone in middle */
  nfr= numberof(fraction);
  rorzc= ror(zcen);

  /* save time */
  times(irec) = time;

  /* get choices for ablation temperatures */
  abltemp= fraction*s(mid,lmax);
  epc= e(mid,2:lmax-1);

  i= prevAbl-1;  /* epc indices 1 less than e indices */
  for (j=nfr ; j>0 ; j--) {
    ablt= abltemp(j);
    while (epc(i)>ablt) i--;  /* find first zone with lower temperature */
    if (j==nfr) prevAbl= i+1;
    rorab(j,irec)= rorzc(i) +
      (ablt-epc(i))*(rorzc(i+1)-rorzc(i))/(epc(i+1)-epc(i));
  }


  /* get peak density and its location */
  dmax(irec)= d(mid, max);

  rordx(irec)= peak(rorzc, d(mid,2:));


  /* get location of peak phot heating */

  rorhx(irec) = peak(rorzc(1:35), h(mid,2:36))


  /* record position, velocity, and acceleration of zone at peak density */
  i= d(mid, mxx);
  zs(irec)= 0.5*(zt(mid,i)+zt(mid,i-1));
  vel(irec)= 0.5*(v(mid,i)+v(mid,i-1));
  acc(irec)= (0.5*(p(mid,i+1)+p(mid,i-1))-p(mid,i))/(ror(i)-ror(i-1));


  /* get density scale height, location of steepest gradient */
  i= min(36, prevScale+4);

  rscl= deriv(d(mid,2:i), rorzc(1:i-1));
  scl(irec)= 1.0/(max(rscl)+1.e-35);

  rorsx(irec) = peak(ror(2:i-1), rscl);
  prevScale= rscl(mxx)+1;

  return;
}

/* Mock up realistic ablation problem.  */

ntimes= 84;
kmax= 11;
lmax= 38;  /* this cannot be changed */

time= 1.5;

fraction= [.4,.5,.6,.7,.8];  /* MUST be in increasing order */

/* the following values represent a very rough fit to an "actual"
   ablating slab problem (that is, the overall shapes of the curves
   are realistic, although the numerical values are not consistent) */
x= span(1,38,38);

y= 4.e-4*atan(x/8) + 1.7e-3*atan((x-16)/3);

y= (2./pi)*(y - y(1));    /* rho*r coordinate */
ror= y(37)-y;             /* ror measured from right (ablating) side */

/* v and p were originally fit on a uniform rho*r grid */

yeven= span(min(y), max(y), 38);
tmp= array(0.0, 38);
tmp(1:35)= -0.2 + 0.4*(x(1:35)-x(1))/(x(35)-x(1));
tmp(35:38)= 0.2 + 0.3*(x(35:38)-x(35))/(x(38)-x(35));
v= interp(tmp, yeven, y)(-:1:kmax,);

/* original fits to density, pressure, temperature, and heating were for
   point centered data -- the strict uncp algorithm will not work here,
   however, since the input was not actually formed as pairwise averages
   of another array... */

func uncp38 (xpc)
{
  x= array(0.0, kmax, 38);
  /* here is "true" uncp algorithm:
     x(2,2)= xpc(1);
     for (i=3 ; i<=38 ; i++) x(2,i)= 2*xpc(i-1)-x(2,i-1);
     x(3:,)= x(2,-,);
   */
  x(2:,2:)= xpc(-,zcen);
  return x;
}


dpc= exp(-40./(x+1)^2 - 50./(x-41)^2 + 1.0-0.22*x);
d= uncp38(dpc);

tmp= -0.01/(x+7)^2 - 0.0001/(x-41)^2;
tmp= tmp-tmp(1) + (tmp(1)-tmp(38))*(x-x(1))/(x(38)-x(1));
tmp(2:10)= tmp(2) + (tmp(10)-tmp(2))*(x(2:10)-x(2))/(x(10)-x(2));
ppc= interp(tmp, yeven, y);
p= uncp38(ppc);

tmp= ppc/dpc;
tmp(19:38)= tmp(19)*(1.0 + 0.1*(x(19:38)-x(19))/(x(38)-x(19)));

tmp(1:5)= max(tmp)/25. + (tmp(5)-max(tmp)/25.)*(x(1:5)-x(1))/(x(5)-x(1));
e= s= uncp38(tmp);

tmp(1:9)= 3.6*(x(1:9)-x(1))/(x(9)-x(1));
tmp(9:20)= 3.6 - 2.6*(x(9:20)-x(9))/(x(20)-x(9));
tmp(21:33)= tmp(20);
tmp(33:38)= 1.0 - (x(33:38)-x(33))/(x(38)-x(33));
h= uncp38(tmp);

/* declare arrays which will hold results of the calculation */
times= array(0.0, ntimes);

rorab= array(0.0, numberof(fraction), ntimes);
rordx= rorhx= rorsx= acc= vel= scl= zs= dmax= 0.0*times;

/* work backwards from d and rhor*r to get zt coordinate */
zt= array(0.0, kmax, 38);
zt(,2:)= (y(dif)/d(2,2:))(-,psum);

/* compute parabolic maximum of a curve

     peak = value of xin where yin is maximum   */

func peak (xin,yin)
{
  imax= yin(mxx);

  if (imax==1 || imax==numberof(yin)) return double(imax);
  ymax= yin(imax);
  dy1= ymax-yin(imax-1);
  dy2= yin(imax+1)-ymax;
  xmax= xin(imax);
  dx1= xmax-xin(imax-1);
  dx2= xin(imax+1)-xmax;
  return xmax - (dy2*dx1*dx1+dy1*dx2*dx2)/(dy2*dx1-dy1*dx2);
}


/*   deriv(y,x) = derivative of y wrt x
  example:

    acc = deriv( deriv(position,time) , time)

       puts acc equal to the second derivative of position wrt time  */

func deriv (y,x)
{
  return y(dif)/x(dif);
}