/*
* $Id: trans.c,v 1.1.1.1 2005/09/18 22:04:57 dhmunro Exp $
* Routines for transporting a ray through a cylindrical mesh.
*/
/* Copyright (c) 2005, The Regents of the University of California.
* All rights reserved.
* This file is part of yorick (http://yorick.sourceforge.net).
* Read the accompanying LICENSE file for details.
*/
#include "trans.h"
/* includes track.h and bound.h as well */
extern double exp(double);
extern double sqrt(double);
extern double fabs(double);
#define SQ(x) ((x)*(x))
#define FUZZ 1.0e-99
#ifdef NLTSS
/* Let Hybrid-C compiler find vector versions of exp, sqrt, etc. */
#include <lcc.h>
#endif
/* Also a CYCRED switch below, originally intended for Crays, but
optimized incorrectly there when last checked. */
/* 1-exp(-x) != x for small x. For x<ONE_MINUS_EXP, just use x. */
#define ONE_MINUS_EXP 1.0e-4
/* ---------------------------------------------------------------------- */
/* Solve the transport equation assuming that both the opacity and
* source function are constant within each zone. Return the overall
* absorption and self-emission along the ray.
*/
void
FlatSource(double *opac, double *source, /* opacity, emissivity */
long kxlm, long ngroup, /* opac, source dimens */
RayPath *path, /* ray path through mesh */
double *absorb, double *selfem, /* attenuation factor and
* self-emission for
* entire path */
double *work) /* scratch space */
{
/* Formulas:
* The transport equation is d(inu)/ds= (knu)(jnu - inu), or
* d(inu)/d(tau)= jnu - inu.
* In a region where knu and jnu are constant, the solution is
* inu= exp(-tau)*inu0 + (1-exp(-tau))*jnu.
*
* Thus, integrating inu along the path requires solving the
* first order recurrence:
* x[i+1]= a[i]*x[i] + b[i]
* where x is inu, a is exp(-tau), and b is (1-exp(-tau))*jnu.
*
* The solution to this recurrence is:
* inu= (absorb)*inu0 + (selfem)
* where absorb is the product of all a[i], and selfem is x[last]
* with x[0]=0. This total attenuation factor and self-emission
* specific intensity (in same units as input source function) are
* the outputs of FlatSource.
*/
long nzones= path->ncuts-1;
long *zone= path->zone;
double *ds= path->ds;
double *tau= work;
double *atten= tau+nzones;
double *jnu= atten+nzones;
long n,i;
register long j;
if (nzones<1) {
if (absorb && selfem) {
for (n=0 ; n<ngroup ; n++) {
*absorb++= 1.0;
*selfem++= 0.0;
}
}
return;
}
/* Eventually, may want to branch to a separate algorithm if the
* number of groups is very large, but the number of edge crossings
* is small...
* For now, the slightly less efficient inner loop on crossings
* is a better choice, since the number of groups is often quite
* small, while the number of crossings is rarely small.
*/
for (n=0 ; n<ngroup ; n++) {
/* Gather the opacity and source function along the ray path.
* Assume that opac[0]==0 and source[0]==0, so no special coding
* required for non-contiguous sections of the path.
*/
/* #pragma nohazard */
for (i=0 ; i<nzones ; i++) {
j= zone[i];
tau[i]= opac[j]*ds[i];
atten[i]= exp(-tau[i]);
jnu[i]= source[j];
}
/* Replace jnu by appropriate "b" function.
* Note that (1-exp(-tau))*jnu --> tau*jnu as tau-->0 may be
* important even for very small tau... */
/* #pragma nohazard */
for (i=0 ; i<nzones ; i++) {
if (fabs(tau[i]) > ONE_MINUS_EXP) jnu[i]*= (1.0-atten[i]);
else jnu[i]*= tau[i];
}
/* Solve the first order recurrence. */
Reduce(atten,jnu,nzones); /* clobbers atten, jnu arrays... */
*absorb= atten[0];
*selfem= jnu[0];
absorb++;
selfem++;
opac+= kxlm;
source+= kxlm;
}
}
/* ---------------------------------------------------------------------- */
/* Solve the transport equation assuming that the source function
* varies linearly with optical depth across a zone. The opacity is
* assumed to be constant within each zone. (Hence, opac must be zone
* centered, while source is point centered.) Return the overall
* absorption and self-emission along the ray.
*/
void
LinearSource(double *opac, double *source, /* opacity, emissivity */
long kxlm, long ngroup, /* opac, source dimens */
RayPath *path, /* ray path through mesh */
double *absorb, double *selfem, /* attenuation factor and
* self-emission for
* entire path */
double *work) /* scratch space */
{
/* Formulas:
* The transport equation is d(inu)/ds= (knu)(jnu - inu), or
* d(inu)/d(tau)= jnu - inu.
* In a region where knu is constant and jnu linear, the solution is
* inu1= exp(-tau)*inu0 + (1-exp(-tau))*jnu0 +
* (1-(1-exp(-tau))/tau)*(jnu1-jnu0).
*
* Thus, integrating inu along the path requires solving the
* first order recurrence:
* x[i+1]= a[i]*x[i] + b[i]
* where x is inu, a is exp(-tau), and b is (1-exp(-tau))*jnu0+
* (1-(1-exp(-tau))/tau)*(jnu1-jnu0).
*
* The solution to this recurrence is:
* inu= (absorb)*inu0 + (selfem)
* where absorb is the product of all a[i], and selfem is x[last]
* with x[0]=0. This total attenuation factor and self-emission
* specific intensity (in same units as input source function) are
* the outputs of LinearSource.
*/
long nzones= path->ncuts-1;
long *zone= path->zone;
double *ds= path->ds;
long *pt1= path->pt1;
long *pt2= path->pt2;
double *f= path->f;
double fi= path->fi;
double ff= path->ff;
double *tau= work;
double *atten= tau+nzones;
double *jnu= atten+nzones; /* jnu has nzones+1 elements! */
long n,i;
register long j,k;
register double f1,f2;
if (nzones<1) {
if (absorb && selfem) {
for (n=0 ; n<ngroup ; n++) {
*absorb++= 1.0;
*selfem++= 0.0;
}
}
return;
}
/* Eventually, may want to branch to a separate algorithm if the
* number of groups is very large, but the number of edge crossings
* is small...
* For now, the slightly less efficient inner loop on crossings
* is a better choice, since the number of groups is often quite
* small, while the number of crossings is rarely small.
*/
for (n=0 ; n<ngroup ; n++) {
/* Gather the opacity and source function along the ray path.
* Assume that opac[0]==0 and source[0]==0, so no special coding
* required for non-contiguous sections of the path.
*/
/* #pragma nohazard */
for (i=0 ; i<nzones ; i++) {
j= zone[i];
tau[i]= opac[j]*ds[i];
atten[i]= exp(-tau[i]);
}
/* Interpolate source to crossing point on each edge */
/* #pragma nohazard */
for (i=0 ; i<=nzones ; i++) {
j= pt1[i]; k= pt2[i];
jnu[i]= (0.5-f[i])*source[j] + (0.5+f[i])*source[k];
}
/* Interpolate source into interior of first and last zone */
f1= (1.0-fi)*jnu[0]+fi*jnu[1]; /* nzones=1 is possible... */
jnu[nzones]= ff*jnu[nzones-1]+(1.0-ff)*jnu[nzones];
jnu[0]= f1;
/* Replace jnu by appropriate "b" function.
* Note that (1-exp(-tau)-tau)*jnu --> (tau^2/2)*jnu as tau-->0
* may be important even for very small tau...
* Vectorizing algorithm requires small positive random number
* in order to make zero divide very unlikely. Magnitude is
* chosen in order to have sub-roundoff effect on
* tau > ONE_MINUS_EXP. */
#define RANDOM_SMALL_TAU (1.5261614e-20*ONE_MINUS_EXP)
/* Note: this vectorizes on Cray as of 3/jul/90 */
/* #pragma nohazard */
for (i=0 ; i<nzones ; i++) {
f1= -atten[i]+(1-atten[i])/(tau[i]+RANDOM_SMALL_TAU);
f2= 1.0 -(1-atten[i])/(tau[i]+RANDOM_SMALL_TAU);
if (fabs(tau[i]) > ONE_MINUS_EXP) {
jnu[i]= f1*jnu[i]+f2*jnu[i+1];
} else {
jnu[i]= 0.5*tau[i]*(jnu[i]+jnu[i+1]);
}
}
/* Solve the first order recurrence. */
Reduce(atten,jnu,nzones); /* clobbers atten, jnu arrays... */
*absorb= atten[0];
*selfem= jnu[0];
absorb++;
selfem++;
opac+= kxlm;
source+= kxlm;
}
}
/* ---------------------------------------------------------------------- */
/* Reduce x[i+1]= a[i]*x[i]+b[i] for i=0,...,n-1 to
* x[n]= a[0]*x[0]+b[0] -i.e.- return a[0], b[0]
* May use an algorithm which clobbers a[1:n-1], b[1:n-1] as well.
*/
void Reduce(double *a, double *b, long n)
{
long i;
#ifndef CYCRED
/* on a non-vector machine, use the straightforward algorithm */
register double tmpa,tmpb;
tmpa= a[0]; tmpb= b[0]; /* miscoded tmpa*b[0] originally?? */
for (i=1 ; i<n ; i++) {
tmpa*= a[i];
tmpb= tmpb*a[i]+b[i];
}
a[0]= tmpa; b[0]= tmpb;
#else
/* on vector machine, a cyclic reduction algorithm is faster */
while (n > 2) { /* log2(n) reduction passes */
/* #pragma nohazard */
for (i=0 ; i<n/2 ; i++) {
a[i]= a[2*i]*a[2*i+1];
b[i]= b[2*i]*a[2*i+1]+b[2*i+1];
}
if (n & 01) { /* keep final odd point if present */
a[i]= a[2*i];
b[i]= b[2*i];
n= i+1;
} else n= i;
}
if (n == 2) {
a[0]= a[0]*a[1];
b[0]= b[0]*a[1]+b[1];
}
#endif
}
/* ---------------------------------------------------------------------- */
/* The LinearSource transport solver requires a point centered source
* function; the given source function is zone centered.
* This point centering scheme is local to the four zones surrounding
* the point, with weights which are large for zones whose optical
* depth (along either k or l) is near 1.
* The algorithm requires 4*(kmax*(lmax+1)+1) doubles as working space.
*
* An unrelated Milne condition is applied for points which lie on
* the vacuum boundary of the problem.
*/
void PtCenterSource(double *opac, double *source, long kxlm, long ngroup,
Mesh *mesh, Boundary *boundary, double *work)
/* Historical note: This algorithm should be algebraically equivalent
* to the algorithm used in DIRT or TDG, excepting a case in the
* boundary correction involving a boundary between non-void zones,
* which can, in fact, never occur.
*/
{
long klmax= mesh->klmax;
long kmax= mesh->kmax;
double *z= mesh->z;
double *r= mesh->r;
int *ireg= mesh->ireg;
/* partition the working space */
double *wp= work; /* wp only needs kmax*lmax... */
double *wk= wp+klmax+kmax+1;
double *wl= wk+klmax+kmax+1;
double *b= wl+klmax+kmax+1;
long i,n;
register double area,dddk,dddl;
register double w00,w01,w10,w11;
/* applying of Milne condition at boundaries logically complex... */
long npoints= boundary->npoints;
long *zone= boundary->zone;
int *side= boundary->side;
/* zone+pt1[side] --> zone+pt2[side] is boundary edge */
long pt1index[4]= { -1L, -1L /* -kmax */, 0L /* -kmax */, 0L };
long pt2index[4]= { 0L, -1L, -1L /* -kmax */, 0L /* -kmax */ };
/* zone+zon[side] is zone interior to boundary by one more than zone */
long zonindex[4]= { 0L /* -kmax */, +1L, 0L /* +kmax */, -1L };
double b1,w1,milne1,b2,w2,milne2,dx,dtau,bx,dtaux,b1b,w1b,milne1b;
long j1,j2,jz,jx,j1b;
int beginSection;
/* without this, not obvious that these are initialized */
b1= w1= milne1= b2= w2= milne2= b1b= w1b= milne1b= 0.0;
j2= j1b= 0;
pt1index[1]-= kmax;
pt1index[2]-= kmax;
pt2index[2]-= kmax;
pt2index[3]-= kmax;
zonindex[0]-= kmax;
zonindex[2]+= kmax;
/* zero workspace in guard zones */
for (i=0 ; i<=kmax ; i++) wk[i]= wl[i]= 0.0;
for (i=klmax ; i<klmax+kmax+1 ; i++) wk[i]= wl[i]= b[i]= 0.0;
for (n=0 ; n<ngroup ; n++) {
/* Compute weighting factors -- zone centered, but separate
* values wk, wl for the delta-k and delta-l directions. */
/* #pragma nohazard */
for (i=kmax+1 ; i<klmax ; i++) {
/* area is actually twice zonal area */
area= (z[i]-z[i-kmax-1])*(r[i-1]-r[i-kmax]) -
(z[i-1]-z[i-kmax])*(r[i]-r[i-kmax-1]);
/* dddk, dddl actually twice median lengths */
dddk= sqrt( SQ(z[i]-z[i-1] + z[i-kmax]-z[i-kmax-1]) +
SQ(r[i]-r[i-1] + r[i-kmax]-r[i-kmax-1]));
dddl= sqrt( SQ(z[i]-z[i-kmax] + z[i-1]-z[i-kmax-1]) +
SQ(r[i]-r[i-kmax] + r[i-1]-r[i-kmax-1]));
if (ireg[i]) {
wp[i]= 1.0/(opac[i]*area+FUZZ);
wk[i]= wp[i]*SQ(1.0-exp(-0.5*opac[i]*dddk));
wl[i]= wp[i]*SQ(1.0-exp(-0.5*opac[i]*dddl));
} else {
wp[i]= wk[i]= wl[i]= 0.0;
}
}
/* Copy source to temporary to add guard zones */
/* #pragma nohazard */
for (i=0 ; i<klmax ; i++) b[i]= source[i];
/* Point centered source is weighted average over
* the four surrounding edges, each of which is centered
* between its two bounding zones. */
/* #pragma nohazard */
for (i=0 ; i<klmax ; i++) {
w00= wk[i]+ wl[i];
w10= wk[i+1]+ wl[i+1];
w01= wk[i+kmax]+ wl[i+kmax];
w11= wk[i+kmax+1]+wl[i+kmax+1];
b[i]= (b[i]*w00+b[i+1]*w10+b[i+kmax]*w01+b[i+kmax+1]*w11) /
(w00+w01+w10+w11+FUZZ);
}
/* Instead of interior weighting, use Milne condition at
* (vacuum) boundaries.
* Note that r=0 or z=0 symmetry boundaries are not in the list,
* because FindBoundaryPoints will not find them. TrimBoundary
* must be used to eliminate khold, lhold, and norad boundaries. */
/* Note that this loop is usually much shorter than the others. */
beginSection= 1;
for (i=0 ; i<npoints ; i++) {
/* Get point indices j1, j2, zone index jz */
jz= zone[i];
if (jz) {
if (beginSection) {
j1= jz + pt1index[side[i]];
} else {
j1= j2;
b1= b2; w1= w2; milne1= milne2;
}
j2= jz + pt2index[side[i]];
b2= source[jz];
if (side[i]&01) w2= wk[jz];
else w2= wl[jz];
/* dtau is half the optical depth of the zone perpendicular
* to the boundary edge -- i.e.- the optical depth from
* zone center to boundary. The Milne factor reduces the
* source so a linear extrapolation reaches 0.0 at
* optical depth 2/3 beyond the zone boundary. If dtau>1
* and the source function is decreasing toward the boundary,
* then extrapolate the source function to optical depth 1
* before applying the Milne reduction. */
dx= sqrt(SQ(z[j2]-z[j1])+SQ(r[j2]-r[j1]));
dtau= 0.25/(dx*wp[jz]+FUZZ);
if (dtau>1.0 && ireg[jx= jz+zonindex[side[i]]]) {
bx= source[jx];
dtaux= 0.25/(dx*wp[jx]+FUZZ);
if (bx > b2) { /* flux (gradient) out of boundary */
if (b2*(dtau+dtaux) > (bx-b2)*(dtau+(2./3.)))
/* effectively extrapolate to optical depth 1 first */
dtau= 1.0+(5./3.)*(bx-b2)*(dtau-1.0)/
(b2*(dtau+dtaux)-(bx-b2)*(dtau-1.0));
} else dtau= 1.0; /* Huh? At least this is continuous... */
}
milne2= 1.0+1.5*dtau;
/* Source at boundary point is weighted mean between the
* two adjacent boundary edges. */
if (beginSection) {
beginSection= 0;
j1b= j1;
b1b= b2; w1b= w2; milne1b= milne2;
b[j1]= b2 / (milne2+FUZZ); /* if closed boundary, changed later */
} else {
b[j1]= (b1*w1+b2*w2) / (w1*milne1+w2*milne2+FUZZ);
}
} else if (!beginSection) { /* final point of boundary section */
beginSection= 1;
if (j2 == j1b) { /* handle closed boundary */
b[j2]= (b1b*w1b+b2*w2) / (w1b*milne1b+w2*milne2+FUZZ);
} else { /* open boundary */
b[j2]= b2 / (milne2+FUZZ);
}
}
}
/* Copy point centered source over zone centered input */
/* #pragma nohazard */
for (i=0 ; i<klmax ; i++) source[i]= b[i];
source+= kxlm;
opac+= kxlm;
}
}
/* ---------------------------------------------------------------------- */
syntax highlighted by Code2HTML, v. 0.9.1