[cig-commits] r18798 - mc/3D/CitcomS/branches/eheien_dev/lib
emheien at geodynamics.org
emheien at geodynamics.org
Fri Jul 22 14:59:35 PDT 2011
Author: emheien
Date: 2011-07-22 14:59:35 -0700 (Fri, 22 Jul 2011)
New Revision: 18798
Modified:
mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c
mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c
mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c
mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c
mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h
mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h
Log:
More work on simplifying and improving tracer code
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c 2011-07-22 20:35:08 UTC (rev 18797)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Composition_related.c 2011-07-22 21:59:35 UTC (rev 18798)
@@ -318,7 +318,7 @@
flavor = i + 1;
E->composition.comp_el[j][i][e] =
E->trace.ntracer_flavor[j][flavor][e] / (double)numtracers;
- fprintf(E->trace.fpt, "%d %d %d %f\n", j, i, e, E->composition.comp_el[j][i][e]);
+ //fprintf(E->trace.fpt, "%d %d %d %f\n", j, i, e, E->composition.comp_el[j][i][e]);
}
}
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c 2011-07-22 20:35:08 UTC (rev 18797)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Full_tracer_advection.c 2011-07-22 21:59:35 UTC (rev 18798)
@@ -74,10 +74,9 @@
double x0, double y0, double z0);
static void crossit(double *cross, double *A, double *B);
static void fix_radius(struct All_variables *E,
- double *radius, double *theta, double *phi,
- double *x, double *y, double *z);
+ SphericalCoord &sc,
+ CartesianCoord &cc);
static void fix_angle(double *angle);
-static void fix_theta_phi(double *theta, double *phi);
static int iget_radial_element(struct All_variables *E,
int j, int iel,
double rad);
@@ -90,7 +89,7 @@
int isend[13][13], double *send[13][13]);
void pdebug(struct All_variables *E, int i);
int full_icheck_cap(struct All_variables *E, int icap,
- double x, double y, double z, double rad);
+ CartesianCoord cc, double rad);
@@ -693,23 +692,22 @@
for (kk=0;kk<irec[j];kk++) {
Tracer new_tracer;
+ CartesianCoord cc;
+ SphericalCoord sc;
ireceive_position=kk*new_tracer.size();
new_tracer.readFromMem(&REC[j][ireceive_position]);
- theta=new_tracer.theta();
- phi=new_tracer.phi();
- rad=new_tracer.rad();
- x=new_tracer.x();
- y=new_tracer.y();
- z=new_tracer.z();
+ sc = new_tracer.getSphericalPos();
+ cc = new_tracer.getCartesianPos();
+
+ iel=(E->trace.iget_element)(E,j,-99,cc,sc);
- iel=(E->trace.iget_element)(E,j,-99,x,y,z,theta,phi,rad);
-
if (iel<1) {
fprintf(E->trace.fpt,"Error(lost souls) - element not here?\n");
- fprintf(E->trace.fpt,"x,y,z-theta,phi,rad: %f %f %f - %f %f %f\n",x,y,z,theta,phi,rad);
+ fprintf(E->trace.fpt,"x,y,z-theta,phi,rad: %f %f %f - %f %f %f\n",
+ cc._x,cc._y,cc._z,sc._theta,sc._phi,sc._rad);
fflush(E->trace.fpt);
exit(10);
}
@@ -750,23 +748,19 @@
int ithatcap, icheck;
int isend_position, ipos;
int lev = E->mesh.levmax;
- double theta, phi, rad;
- double x, y, z;
TracerList::iterator tr;
+ CartesianCoord cc;
/* transfer tracers from rlater to send */
for (tr=E->trace.escaped_tracers[j].begin();tr!=E->trace.escaped_tracers[j].end();++tr) {
- rad=tr->rad();
- x=tr->x();
- y=tr->y();
- z=tr->z();
+ cc = tr->getCartesianPos();
/* first check same cap if nprocz>1 */
if (E->parallel.nprocz>1) {
ithatcap=0;
- icheck=full_icheck_cap(E,ithatcap,x,y,z,rad);
+ icheck=full_icheck_cap(E,ithatcap,cc,tr->rad());
if (icheck==1) goto foundit;
}
@@ -775,7 +769,7 @@
for (pp=1;pp<=E->parallel.TNUM_PASS[lev][j];pp++) {
ithatcap=pp;
- icheck=full_icheck_cap(E,ithatcap,x,y,z,rad);
+ icheck=full_icheck_cap(E,ithatcap,cc,tr->rad());
if (icheck==1) goto foundit;
}
@@ -783,8 +777,8 @@
/* should not be here */
if (icheck!=1) {
fprintf(E->trace.fpt,"Error(lost souls)-should not be here\n");
- fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",x,y,z,rad);
- icheck=full_icheck_cap(E,0,x,y,z,rad);
+ fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",cc._x,cc._y,cc._z,tr->rad());
+ icheck=full_icheck_cap(E,0,cc,tr->rad());
if (icheck==1) fprintf(E->trace.fpt," icheck here!\n");
else fprintf(E->trace.fpt,"icheck not here!\n");
fflush(E->trace.fpt);
@@ -857,9 +851,6 @@
const double eps=-1e-4;
- void sphere_to_cart();
-
-
/* find u and v using spherical coordinates */
spherical_to_uv(E,j,theta,phi,&u,&v);
@@ -909,7 +900,7 @@
sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
ival=icheck_element(E,j,nelem,x,y,z,rad);
fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
- ival=(E->trace.iget_element)(E,j,-99,x,y,z,theta,phi,rad);
+ ival=(E->trace.iget_element)(E,j,-99,CartesianCoord(x,y,z),SphericalCoord(theta,phi,rad));
fprintf(E->trace.fpt,"New Element?: %d\n",ival);
ival=icheck_column_neighbors(E,j,nelem,x,y,z,rad);
fprintf(E->trace.fpt,"New Element (neighs)?: %d\n",ival);
@@ -994,10 +985,9 @@
/* 5 6 5 7 */
/* 6 7 6 8 */
-void full_get_velocity(struct All_variables *E,
+CartesianCoord full_get_velocity(struct All_variables *E,
int j, int nelem,
- double theta, double phi, double rad,
- double *velocity_vector)
+ SphericalCoord sc)
{
int iwedge;
const int sphere_key = 0;
@@ -1006,9 +996,7 @@
double VV[4][9];
double vx[7],vy[7],vz[7];
- void velo_from_element_d();
-
- full_get_shape_functions(E, shape, nelem, theta, phi, rad);
+ full_get_shape_functions(E, shape, nelem, sc._theta, sc._phi, sc._rad);
iwedge=shape[0];
/* get cartesian velocity */
@@ -1059,16 +1047,12 @@
vz[6]=VV[3][8];
}
- velocity_vector[1]=vx[1]*shape[1]+vx[2]*shape[2]+shape[3]*vx[3]+
- vx[4]*shape[4]+vx[5]*shape[5]+shape[6]*vx[6];
- velocity_vector[2]=vy[1]*shape[1]+vy[2]*shape[2]+shape[3]*vy[3]+
- vy[4]*shape[4]+vy[5]*shape[5]+shape[6]*vy[6];
- velocity_vector[3]=vz[1]*shape[1]+vz[2]*shape[2]+shape[3]*vz[3]+
- vz[4]*shape[4]+vz[5]*shape[5]+shape[6]*vz[6];
-
-
-
- return;
+ return CartesianCoord(vx[1]*shape[1]+vx[2]*shape[2]+shape[3]*vx[3]+
+ vx[4]*shape[4]+vx[5]*shape[5]+shape[6]*vx[6],
+ vy[1]*shape[1]+vy[2]*shape[2]+shape[3]*vy[3]+
+ vy[4]*shape[4]+vy[5]*shape[5]+shape[6]*vy[6],
+ vz[1]*shape[1]+vz[2]*shape[2]+shape[3]*vz[3]+
+ vz[4]*shape[4]+vz[5]*shape[5]+shape[6]*vz[6]);
}
/***************************************************************/
@@ -2199,7 +2183,7 @@
/* a given cap */
/* */
int full_icheck_cap(struct All_variables *E, int icap,
- double x, double y, double z, double rad)
+ CartesianCoord cc, double rad)
{
double test_point[4];
@@ -2228,9 +2212,9 @@
/* test_point - project to outer radius */
- test_point[1]=x/rad;
- test_point[2]=y/rad;
- test_point[3]=z/rad;
+ test_point[1]=cc._x/rad;
+ test_point[2]=cc._y/rad;
+ test_point[3]=cc._z/rad;
ival=icheck_bounds(E,test_point,rnode[1],rnode[2],rnode[3],rnode[4]);
@@ -2432,39 +2416,36 @@
/* This function moves particles back in bounds if they left */
/* during advection */
-static void fix_radius(struct All_variables *E,
- double *radius, double *theta, double *phi,
- double *x, double *y, double *z)
+static void fix_radius(struct All_variables *E, SphericalCoord &sc, CartesianCoord &cc)
{
double sint,cost,sinf,cosf,rad;
double max_radius, min_radius;
-
+ int changed = 0;
+
max_radius = E->sphere.ro - E->trace.box_cushion;
min_radius = E->sphere.ri + E->trace.box_cushion;
-
- if (*radius > max_radius) {
- *radius=max_radius;
+
+ if (sc._rad > max_radius) {
+ changed = 1;
+ sc._rad=max_radius;
rad=max_radius;
- cost=cos(*theta);
- sint=sqrt(1.0-cost*cost);
- cosf=cos(*phi);
- sinf=sin(*phi);
- *x=rad*sint*cosf;
- *y=rad*sint*sinf;
- *z=rad*cost;
}
- if (*radius < min_radius) {
- *radius=min_radius;
+ if (sc._rad < min_radius) {
+ changed = 1;
+ sc._rad=min_radius;
rad=min_radius;
- cost=cos(*theta);
- sint=sqrt(1.0-cost*cost);
- cosf=cos(*phi);
- sinf=sin(*phi);
- *x=rad*sint*cosf;
- *y=rad*sint*sinf;
- *z=rad*cost;
}
-
+
+ if (changed) {
+ cost=cos(sc._theta);
+ sint=sqrt(1.0-cost*cost);
+ cosf=cos(sc._phi);
+ sinf=sin(sc._phi);
+ cc._x=rad*sint*cosf;
+ cc._y=rad*sint*sinf;
+ cc._z=rad*cost;
+ }
+
return;
}
@@ -2487,30 +2468,6 @@
return;
}
-/******************************************************************/
-/* FIX THETA PHI */
-/* */
-/* This function constrains the value of theta to be */
-/* between 0 and PI, and */
-/* this function constrains the value of phi to be */
-/* between 0 and 2 PI */
-/* */
-static void fix_theta_phi(double *theta, double *phi)
-{
- const double two_pi=2.0*M_PI;
-
- fix_angle(theta);
-
- if (*theta > M_PI) {
- *theta = two_pi - *theta;
- *phi += M_PI;
- }
-
- fix_angle(phi);
-
- return;
-}
-
/********** IGET ELEMENT *****************************************/
/* */
/* This function returns the the real element for a given point. */
@@ -2521,8 +2478,8 @@
int full_iget_element(struct All_variables *E,
int j, int iprevious_element,
- double x, double y, double z,
- double theta, double phi, double rad)
+ CartesianCoord cc,
+ SphericalCoord sc)
{
int icheck_processor_shell();
int iregel;
@@ -2541,14 +2498,13 @@
ely=E->lmesh.ely;
elz=E->lmesh.elz;
-
ntheta=0;
nphi=0;
/* check the radial range */
if (E->parallel.nprocz>1)
{
- ival=icheck_processor_shell(E,j,rad);
+ ival=icheck_processor_shell(E,j,sc._rad);
if (ival!=1) return -99;
}
@@ -2560,7 +2516,7 @@
/* get regular element number */
- iregel=iget_regel(E,j,theta,phi,&ntheta,&nphi);
+ iregel=iget_regel(E,j,sc._theta,sc._phi,&ntheta,&nphi);
if (iregel<=0)
{
return -99;
@@ -2579,7 +2535,7 @@
if (iprevious_element>0)
{
- ival=icheck_element_column(E,j,iprevious_element,x,y,z,rad);
+ ival=icheck_element_column(E,j,iprevious_element,cc._x,cc._y,cc._z,sc._rad);
if (ival==1)
{
iel=iprevious_element;
@@ -2600,7 +2556,7 @@
if (nelem!=iprevious_element)
{
- ival=icheck_element_column(E,j,nelem,x,y,z,rad);
+ ival=icheck_element_column(E,j,nelem,cc._x,cc._y,cc._z,sc._rad);
if (ival==1)
{
iel=nelem;
@@ -2617,7 +2573,7 @@
if (iprevious_element>0)
{
- iel=icheck_column_neighbors(E,j,iprevious_element,x,y,z,rad);
+ iel=icheck_column_neighbors(E,j,iprevious_element,cc._x,cc._y,cc._z,sc._rad);
if (iel>0)
{
goto foundit;
@@ -2626,7 +2582,7 @@
/* check if still in cap */
- ival=full_icheck_cap(E,0,x,y,z,rad);
+ ival=full_icheck_cap(E,0,cc,sc._rad);
if (ival==0)
{
return -99;
@@ -2643,7 +2599,7 @@
icorner[4]=elxz*ely;
for (kk=1;kk<=4;kk++)
{
- ival=icheck_element_column(E,j,icorner[kk],x,y,z,rad);
+ ival=icheck_element_column(E,j,icorner[kk],cc._x,cc._y,cc._z,sc._rad);
if (ival>0)
{
iel=icorner[kk];
@@ -2661,7 +2617,7 @@
for (kk=1;kk<=ichoice;kk++)
{
ineighbor=E->trace.regtoel[j][kk][iregel];
- iel=icheck_column_neighbors(E,j,ineighbor,x,y,z,rad);
+ iel=icheck_column_neighbors(E,j,ineighbor,cc._x,cc._y,cc._z,sc._rad);
if (iel>0)
{
goto foundit;
@@ -2675,7 +2631,7 @@
E->trace.istat1++;
- iel=icheck_all_columns(E,j,x,y,z,rad);
+ iel=icheck_all_columns(E,j,cc._x,cc._y,cc._z,sc._rad);
/*
fprintf(E->trace.fpt,"WARNING(full_iget_element)-doing a full search!\n");
@@ -2703,7 +2659,7 @@
fprintf(E->trace.fpt,"Error(full_iget_element) - element not found\n");
fprintf(E->trace.fpt,"x,y,z,theta,phi,iregel %.15e %.15e %.15e %.15e %.15e %d\n",
- x,y,z,theta,phi,iregel);
+ cc._x,cc._y,cc._z,sc._theta,sc._phi,iregel);
fflush(E->trace.fpt);
return -1;
@@ -2711,7 +2667,7 @@
/* find radial element */
- ifinal_iel=iget_radial_element(E,j,iel,rad);
+ ifinal_iel=iget_radial_element(E,j,iel,sc._rad);
return ifinal_iel;
}
@@ -3013,11 +2969,11 @@
/* phi and theta are within the proper degree range. */
void full_keep_within_bounds(struct All_variables *E,
- double *x, double *y, double *z,
- double *theta, double *phi, double *rad)
+ CartesianCoord &cc,
+ SphericalCoord &sc)
{
- fix_theta_phi(theta, phi);
- fix_radius(E,rad,theta,phi,x,y,z);
+ sc.constrainThetaPhi();
+ fix_radius(E,sc,cc);
return;
}
@@ -3307,10 +3263,6 @@
double time;
double path,dtpath;
- void sphere_to_cart();
- void cart_to_sphere();
- void analytical_test_function();
-
theta_0=x0_s[1];
phi_0=x0_s[2];
rad_0=x0_s[3];
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c 2011-07-22 20:35:08 UTC (rev 18797)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Regional_tracer_advection.c 2011-07-22 21:59:35 UTC (rev 18798)
@@ -270,8 +270,8 @@
int regional_iget_element(struct All_variables *E,
int m, int iprevious_element,
- double dummy1, double dummy2, double dummy3,
- double theta, double phi, double rad)
+ CartesianCoord dummy,
+ SphericalCoord sc)
{
int e, i, j, k;
int ii, jj, kk;
@@ -291,9 +291,9 @@
i = (e / elz) % elx;
j = e / (elz*elx);
- ii = isearch_neighbors(E->trace.x_space, elx+1, theta, i);
- jj = isearch_neighbors(E->trace.y_space, ely+1, phi, j);
- kk = isearch_neighbors(E->trace.z_space, elz+1, rad, k);
+ ii = isearch_neighbors(E->trace.x_space, elx+1, sc._theta, i);
+ jj = isearch_neighbors(E->trace.y_space, ely+1, sc._phi, j);
+ kk = isearch_neighbors(E->trace.z_space, elz+1, sc._rad, k);
if (ii>=0 && jj>=0 && kk>=0)
return jj*elx*elz + ii*elz + kk + 1;
@@ -301,9 +301,9 @@
/* Search all elements if either the previous element is unknown */
/* or failed to find in the neighboring elements */
- ii = isearch_all(E->trace.x_space, elx+1, theta);
- jj = isearch_all(E->trace.y_space, ely+1, phi);
- kk = isearch_all(E->trace.z_space, elz+1, rad);
+ ii = isearch_all(E->trace.x_space, elx+1, sc._theta);
+ jj = isearch_all(E->trace.y_space, ely+1, sc._phi);
+ kk = isearch_all(E->trace.z_space, elz+1, sc._rad);
if (ii<0 || jj<0 || kk<0)
return -99;
@@ -368,7 +368,7 @@
/* a given cap */
/* */
int regional_icheck_cap(struct All_variables *E, int icap,
- double theta, double phi, double rad, double junk)
+ SphericalCoord sc, double junk)
{
double theta_min, theta_max;
double phi_min, phi_max;
@@ -382,8 +382,8 @@
phi_min = E->trace.phi_cap[icap][2];
phi_max = E->trace.phi_cap[icap][4];
- if ((theta >= theta_min) && (theta < theta_max) &&
- (phi >= phi_min) && (phi < phi_max))
+ if ((sc._theta >= theta_min) && (sc._theta < theta_max) &&
+ (sc._phi >= phi_min) && (sc._phi < phi_max))
return 1;
//TODO: deal with south west bounds
@@ -480,19 +480,19 @@
/******** GET VELOCITY ***************************************/
-void regional_get_velocity(struct All_variables *E,
+CartesianCoord regional_get_velocity(struct All_variables *E,
int m, int nelem,
- double theta, double phi, double rad,
- double *velocity_vector)
+ SphericalCoord sc)
{
void velo_from_element_d();
double shp[9], VV[4][9], tmp;
int n, d, node;
const int sphere_key = 0;
+ double velocity_vector[4];
/* get shape functions at (theta, phi, rad) */
- regional_get_shape_functions(E, shp, nelem, theta, phi, rad);
+ regional_get_shape_functions(E, shp, nelem, sc._theta, sc._phi, sc._rad);
/* get cartesian velocity */
@@ -530,49 +530,48 @@
fflush(E->trace.fpt);
/**/
- return;
+ return CartesianCoord(velocity_vector[1], velocity_vector[2], velocity_vector[3]);
}
void regional_keep_within_bounds(struct All_variables *E,
- double *x, double *y, double *z,
- double *theta, double *phi, double *rad)
+ CartesianCoord &cc,
+ SphericalCoord &sc)
{
- void sphere_to_cart();
int changed = 0;
- if (*theta > E->control.theta_max - E->trace.box_cushion) {
- *theta = E->control.theta_max - E->trace.box_cushion;
+ if (sc._theta > E->control.theta_max - E->trace.box_cushion) {
+ sc._theta = E->control.theta_max - E->trace.box_cushion;
changed = 1;
}
- if (*theta < E->control.theta_min + E->trace.box_cushion) {
- *theta = E->control.theta_min + E->trace.box_cushion;
+ if (sc._theta < E->control.theta_min + E->trace.box_cushion) {
+ sc._theta = E->control.theta_min + E->trace.box_cushion;
changed = 1;
}
- if (*phi > E->control.fi_max - E->trace.box_cushion) {
- *phi = E->control.fi_max - E->trace.box_cushion;
+ if (sc._phi > E->control.fi_max - E->trace.box_cushion) {
+ sc._phi = E->control.fi_max - E->trace.box_cushion;
changed = 1;
}
- if (*phi < E->control.fi_min + E->trace.box_cushion) {
- *phi = E->control.fi_min + E->trace.box_cushion;
+ if (sc._phi < E->control.fi_min + E->trace.box_cushion) {
+ sc._phi = E->control.fi_min + E->trace.box_cushion;
changed = 1;
}
- if (*rad > E->sphere.ro - E->trace.box_cushion) {
- *rad = E->sphere.ro - E->trace.box_cushion;
+ if (sc._rad > E->sphere.ro - E->trace.box_cushion) {
+ sc._rad = E->sphere.ro - E->trace.box_cushion;
changed = 1;
}
- if (*rad < E->sphere.ri + E->trace.box_cushion) {
- *rad = E->sphere.ri + E->trace.box_cushion;
+ if (sc._rad < E->sphere.ri + E->trace.box_cushion) {
+ sc._rad = E->sphere.ri + E->trace.box_cushion;
changed = 1;
}
if (changed)
- sphere_to_cart(E, *theta, *phi, *rad, x, y, z);
+ cc = sc.toCartesian();
return;
@@ -909,7 +908,7 @@
/* check radius first, since it is cheaper */
inside = icheck_processor_shell(E, j, rad);
if (inside == 1)
- inside = regional_icheck_cap(E, 0, theta, phi, rad, rad);
+ inside = regional_icheck_cap(E, 0, SphericalCoord(theta, phi, rad), rad);
else
inside = 0;
@@ -925,7 +924,7 @@
new_tracer.readFromMem(&recv[ipos]);
/* found the element */
- iel = regional_iget_element(E, j, -99, 0, 0, 0, theta, phi, rad);
+ iel = regional_iget_element(E, j, -99, CartesianCoord(0, 0, 0), SphericalCoord(theta, phi, rad));
if (iel<1) {
fprintf(E->trace.fpt, "Error(regional lost souls) - "
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c 2011-07-22 20:35:08 UTC (rev 18797)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/Tracer_setup.c 2011-07-22 21:59:35 UTC (rev 18798)
@@ -62,9 +62,9 @@
void count_tracers_of_flavors(struct All_variables *E);
int full_icheck_cap(struct All_variables *E, int icap,
- double x, double y, double z, double rad);
+ CartesianCoord cc, double rad);
int regional_icheck_cap(struct All_variables *E, int icap,
- double x, double y, double z, double rad);
+ SphericalCoord sc, double rad);
static void find_tracers(struct All_variables *E);
static void predict_tracers(struct All_variables *E);
@@ -143,6 +143,72 @@
+SphericalCoord CartesianCoord::toSpherical(void) const
+{
+ double temp;
+ double theta, phi, rad;
+
+ temp=_x*_x+_y*_y;
+
+ theta=atan2(sqrt(temp),_z);
+ phi=myatan(_y,_x);
+ rad=sqrt(temp+_z*_z);
+
+ return SphericalCoord(theta, phi, rad);
+}
+
+CartesianCoord SphericalCoord::toCartesian(void) const
+{
+ double sint,cost,sinf,cosf;
+ double x, y, z;
+
+ sint=sin(_theta);
+ cost=cos(_theta);
+ sinf=sin(_phi);
+ cosf=cos(_phi);
+
+ x=_rad*sint*cosf;
+ y=_rad*sint*sinf;
+ z=_rad*cost;
+
+ return CartesianCoord(x,y,z);
+}
+
+
+
+const CartesianCoord CartesianCoord::operator+(const CartesianCoord &other) const {
+ return CartesianCoord( this->_x+other._x,
+ this->_y+other._y,
+ this->_z+other._z);
+}
+
+// Multiply each element by a constant factor
+const CartesianCoord CartesianCoord::operator*(const double &val) const {
+ return CartesianCoord( this->_x*val,
+ this->_y*val,
+ this->_z*val);
+}
+
+
+// Returns a constrained angle between 0 and 2 PI
+double SphericalCoord::constrainAngle(const double angle) const {
+ const double two_pi = 2.0*M_PI;
+
+ return angle - two_pi * floor(angle/two_pi);
+}
+
+// Constrains theta to be between 0 and PI while simultaneously constraining phi between 0 and 2 PI
+void SphericalCoord::constrainThetaPhi(void) {
+ _theta = constrainAngle(_theta);
+ if (_theta > M_PI) {
+ _theta = 2*M_PI - _theta;
+ _phi += M_PI;
+ }
+ _phi = constrainAngle(_phi);
+}
+
+
+
void tracer_input(struct All_variables *E)
{
void full_tracer_input();
@@ -322,8 +388,6 @@
E->trace.advection_time += CPU_time0() - begin_time;
tracer_post_processing(E);
-
- return;
}
@@ -394,160 +458,91 @@
/* */
/* This function predicts tracers performing an euler step */
/* */
-/* */
-/* Note positions used in tracer array */
-/* [positions 0-5 are always fixed with current coordinates */
-/* Positions 6-8 contain original Cartesian coordinates. */
-/* Positions 9-11 contain original Cartesian velocities. */
-/* */
static void predict_tracers(struct All_variables *E)
{
-
- int numtracers;
- int j;
- int nelem;
-
- double dt;
- double theta0,phi0,rad0;
- double x0,y0,z0;
- double theta_pred,phi_pred,rad_pred;
- double x_pred,y_pred,z_pred;
- double velocity_vector[4];
-
+ int j, nelem;
+ double dt;
+ SphericalCoord sc0, sc_pred;
+ CartesianCoord cc0, cc_pred, velocity_vector;
TracerList::iterator tr;
-
-
+
dt=E->advection.timestep;
-
-
+
for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
+
for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
-
- theta0=tr->theta();
- phi0=tr->phi();
- rad0=tr->rad();
- x0=tr->x();
- y0=tr->y();
- z0=tr->z();
-
- nelem=tr->ielement();
- (E->trace.get_velocity)(E,j,nelem,theta0,phi0,rad0,velocity_vector);
-
- x_pred=x0+velocity_vector[1]*dt;
- y_pred=y0+velocity_vector[2]*dt;
- z_pred=z0+velocity_vector[3]*dt;
-
-
- /* keep in box */
-
- cart_to_sphere(E,x_pred,y_pred,z_pred,&theta_pred,&phi_pred,&rad_pred);
- (E->trace.keep_within_bounds)(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
-
- /* Current Coordinates are always kept in positions 0-5. */
-
- tr->setCoords(SphericalCoord(theta_pred, phi_pred, rad_pred), CartesianCoord(x_pred, y_pred, z_pred));
-
- /* Fill in original coords (positions 6-8) */
- /* Fill in original velocities (positions 9-11) */
-
- tr->setOrigVals(CartesianCoord(x0, y0, z0),
- CartesianCoord(velocity_vector[1], velocity_vector[2], velocity_vector[3]));
-
- } /* end kk, predicting tracers */
+ // Get initial tracer position
+ sc0 = tr->getSphericalPos();
+ cc0 = tr->getCartesianPos();
+
+ // Interpolate the velocity of the tracer from element it is in
+ nelem = tr->ielement();
+ velocity_vector = (E->trace.get_velocity)(E,j,nelem,sc0);
+
+ // Advance the tracer location in an Euler step
+ cc_pred = cc0 + velocity_vector * dt;
+
+ // Ensure tracer remains in the boundaries
+ sc_pred = cc_pred.toSpherical();
+ (E->trace.keep_within_bounds)(E, cc_pred, sc_pred);
+
+ // Set new tracer coordinates
+ tr->setCoords(sc_pred, cc_pred);
+
+ // Save original coords and velocities for later correction
+
+ tr->setOrigVals(cc0, velocity_vector);
+
+ } /* end predicting tracers */
} /* end caps */
-
+
/* find new tracer elements and caps */
-
find_tracers(E);
-
- return;
-
}
-
/*********** CORRECT TRACERS **********************************************/
/* */
/* This function corrects tracers using both initial and */
/* predicted velocities */
/* */
-/* */
-/* Note positions used in tracer array */
-/* [positions 0-5 are always fixed with current coordinates */
-/* Positions 6-8 contain original Cartesian coordinates. */
-/* Positions 9-11 contain original Cartesian velocities. */
-/* */
-
static void correct_tracers(struct All_variables *E)
{
-
- int j;
- int nelem;
-
-
- double dt;
- double x0,y0,z0;
- double theta_pred,phi_pred,rad_pred;
- double x_pred,y_pred,z_pred;
- double theta_cor,phi_cor,rad_cor;
- double x_cor,y_cor,z_cor;
- double velocity_vector[4];
- double Vx0,Vy0,Vz0;
- double Vx_pred,Vy_pred,Vz_pred;
-
+ int j, nelem;
+ double dt;
TracerList::iterator tr;
-
-
+ CartesianCoord orig_pos, orig_vel, pred_vel, new_coord;
+ SphericalCoord new_coord_sph;
+
dt=E->advection.timestep;
-
-
+
for (j=1;j<=E->sphere.caps_per_proc;j++) {
for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr) {
-
- theta_pred=tr->theta();
- phi_pred=tr->phi();
- rad_pred=tr->rad();
- x_pred=tr->x();
- y_pred=tr->y();
- z_pred=tr->z();
-
- x0=tr->getOrigCartesianPos()._x;
- y0=tr->getOrigCartesianPos()._y;
- z0=tr->getOrigCartesianPos()._z;
-
- Vx0=tr->getCartesianVel()._x;
- Vy0=tr->getCartesianVel()._y;
- Vz0=tr->getCartesianVel()._z;
-
- nelem=tr->ielement();
-
- (E->trace.get_velocity)(E,j,nelem,theta_pred,phi_pred,rad_pred,velocity_vector);
-
- Vx_pred=velocity_vector[1];
- Vy_pred=velocity_vector[2];
- Vz_pred=velocity_vector[3];
-
- x_cor=x0 + dt * 0.5*(Vx0+Vx_pred);
- y_cor=y0 + dt * 0.5*(Vy0+Vy_pred);
- z_cor=z0 + dt * 0.5*(Vz0+Vz_pred);
-
- cart_to_sphere(E,x_cor,y_cor,z_cor,&theta_cor,&phi_cor,&rad_cor);
- (E->trace.keep_within_bounds)(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
-
- /* Fill in Current Positions (other positions are no longer important) */
-
- tr->setCoords(SphericalCoord(theta_cor, phi_cor, rad_cor), CartesianCoord(x_cor, y_cor, z_cor));
-
- } /* end kk, correcting tracers */
+
+ // Get the original tracer position and velocity
+ orig_pos = tr->getOrigCartesianPos();
+ orig_vel = tr->getCartesianVel();
+
+ // Get the predicted velocity by interpolating in the current element
+ nelem = tr->ielement();
+ pred_vel = (E->trace.get_velocity)(E, j, nelem, tr->getSphericalPos());
+
+ // Correct the tracer position based on the original and new velocities
+ new_coord = orig_pos + (orig_vel+pred_vel) * (dt*0.5);
+ new_coord_sph = new_coord.toSpherical();
+
+ // Ensure tracer remains in boundaries
+ (E->trace.keep_within_bounds)(E, new_coord, new_coord_sph);
+
+ // Fill in current positions (original values are no longer important)
+ tr->setCoords(new_coord_sph, new_coord);
+
+ } /* end correcting tracers */
} /* end caps */
-
+
/* find new tracer elements and caps */
-
find_tracers(E);
-
- return;
}
@@ -560,46 +555,33 @@
static void find_tracers(struct All_variables *E)
{
- int iel;
- int kk;
- int j;
- int iprevious_element;
-
- double theta,phi,rad;
- double x,y,z;
- double time_stat1;
- double time_stat2;
-
+ int iel, j, iprevious_element;
TracerList::iterator tr;
+ double begin_time = CPU_time0();
- double CPU_time0();
- double begin_time = CPU_time0();
-
-
for (j=1;j<=E->sphere.caps_per_proc;j++) {
-
/* initialize arrays and statistical counters */
E->trace.istat1=0;
- for (kk=0;kk<=4;kk++) {
+ for (int kk=0;kk<=4;kk++) {
E->trace.istat_ichoice[j][kk]=0;
}
for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();) {
- theta=tr->theta();
- phi=tr->phi();
- rad=tr->rad();
- x=tr->x();
- y=tr->y();
- z=tr->z();
+ CartesianCoord cc;
+ SphericalCoord sc;
+
+ sc = tr->getSphericalPos();
+ cc = tr->getCartesianPos();
iprevious_element=tr->ielement();
- iel=(E->trace.iget_element)(E,j,iprevious_element,x,y,z,theta,phi,rad);
+ iel=(E->trace.iget_element)(E,j,iprevious_element,cc,sc);
/* debug *
- fprintf(E->trace.fpt,"BB. kk %d %d %d %d %f %f %f %f %f %f\n",kk,j,iprevious_element,iel,x,y,z,theta,phi,rad);
+ fprintf(E->trace.fpt,"BB. kk %d %d %d %f %f %f %f %f %f\n",
+ j,iprevious_element,iel,cc._x,cc._y,cc._z,sc._theta,sc._phi,sc._rad);
fflush(E->trace.fpt);
/**/
@@ -629,11 +611,6 @@
/* Now take care of tracers that exited cap */
- /* REMOVE */
- /*
- parallel_process_termination();
- */
-
if (E->parallel.nprocxy == 12)
full_lost_souls(E);
else
@@ -721,7 +698,6 @@
find_tracers(E);
-
/* count # of tracers of each flavor */
if (E->trace.nflavors > 0)
@@ -756,15 +732,11 @@
generate_random_tracers(E, tracers_cap, j);
-
-
}/* end j */
/* Initialize tracer flavors */
if (E->trace.nflavors) init_tracer_flavors(E);
-
- return;
}
@@ -777,8 +749,6 @@
int number_of_tries=0;
int max_tries;
- double x,y,z;
- double theta,phi,rad;
double xmin,xmax,ymin,ymax,zmin,zmax;
double random1,random2,random3;
@@ -790,6 +760,7 @@
xmin = ymin = zmin = E->sphere.ro;
xmax = ymax = zmax = -E->sphere.ro;
for (kk=1; kk<=E->lmesh.nno; kk++) {
+ double x,y,z;
x = E->x[j][1][kk];
y = E->x[j][2][kk];
z = E->x[j][3][kk];
@@ -826,32 +797,35 @@
random3=(1.0*rand())/(1.0*RAND_MAX);
#endif
- x=xmin+random1*(xmax-xmin);
- y=ymin+random2*(ymax-ymin);
- z=zmin+random3*(zmax-zmin);
+ CartesianCoord new_cc;
+ SphericalCoord new_sc;
+
+ new_cc = CartesianCoord(xmin+random1*(xmax-xmin),
+ ymin+random2*(ymax-ymin),
+ zmin+random3*(zmax-zmin));
/* first check if within shell */
- cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
+ new_sc = new_cc.toSpherical();
- if (rad>=E->sx[j][3][E->lmesh.noz]) continue;
- if (rad<E->sx[j][3][1]) continue;
+ if (new_sc._rad>=E->sx[j][3][E->lmesh.noz]) continue;
+ if (new_sc._rad<E->sx[j][3][1]) continue;
/* check if in current cap */
if (E->parallel.nprocxy==1)
- ival=regional_icheck_cap(E,0,theta,phi,rad,rad);
+ ival=regional_icheck_cap(E,0,new_sc,new_sc._rad);
else
- ival=full_icheck_cap(E,0,x,y,z,rad);
+ ival=full_icheck_cap(E,0,new_cc,new_sc._rad);
if (ival!=1) continue;
/* Made it, so record tracer information */
- (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
+ (E->trace.keep_within_bounds)(E, new_cc, new_sc);
Tracer new_tracer;
- new_tracer.setCoords(SphericalCoord(theta, phi, rad), CartesianCoord(x, y, z));
+ new_tracer.setCoords(new_sc, new_cc);
E->trace.tracers[j].push_back(new_tracer);
} /* end while */
@@ -877,8 +851,6 @@
int icushion;
int i, j;
- double x,y,z;
- double theta,phi,rad;
double buffer[100];
FILE *fptracer;
@@ -900,7 +872,6 @@
exit(10);
}
-
/* initially size tracer arrays to number of tracers divided by processors */
icushion=100;
@@ -915,7 +886,10 @@
allocate_tracer_arrays(E,j,iestimate);
for (kk=1;kk<=number_of_tracers;kk++) {
- int len, ncol;
+ SphericalCoord in_coord_sph;
+ CartesianCoord in_coord_cc;
+ int len, ncol;
+
ncol = 3 + 1; // ERIC: E->trace.number_of_extra_quantities;
len = read_double_vector(fptracer, ncol, buffer);
@@ -925,27 +899,23 @@
exit(10);
}
- theta = buffer[0];
- phi = buffer[1];
- rad = buffer[2];
+ in_coord_sph.readFromMem(buffer);
+ in_coord_cc = in_coord_sph.toCartesian();
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
-
-
/* make sure theta, phi is in range, and radius is within bounds */
- (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
+ (E->trace.keep_within_bounds)(E,in_coord_cc,in_coord_sph);
/* check whether tracer is within processor domain */
icheck=1;
- if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,rad);
+ if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,in_coord_sph._rad);
if (icheck!=1) continue;
if (E->parallel.nprocxy==1)
- icheck=regional_icheck_cap(E,0,theta,phi,rad,rad);
+ icheck=regional_icheck_cap(E,0,in_coord_sph,in_coord_sph._rad);
else
- icheck=full_icheck_cap(E,0,x,y,z,rad);
+ icheck=full_icheck_cap(E,0,in_coord_cc,in_coord_sph._rad);
if (icheck==0) continue;
@@ -953,7 +923,7 @@
Tracer new_tracer;
- new_tracer.setCoords(SphericalCoord(theta, phi, rad), CartesianCoord(x,y,z));
+ new_tracer.setCoords(in_coord_sph, in_coord_cc);
new_tracer.set_flavor(buffer[3]);
E->trace.tracers[j].push_back(new_tracer);
@@ -1018,8 +988,6 @@
FILE *fp1;
-
-
/* deal with different output formats */
#ifdef USE_GZDIR
if(strcmp(E->output.format, "ascii-gz") == 0){
@@ -1069,6 +1037,9 @@
for (kk=1;kk<=numtracers;kk++) {
int len, ncol;
+ SphericalCoord sc;
+ CartesianCoord cc;
+
ncol = 3 + 1;
len = read_double_vector(fp1, ncol, buffer);
@@ -1078,19 +1049,16 @@
exit(10);
}
- theta = buffer[0];
- phi = buffer[1];
- rad = buffer[2];
+ sc.readFromMem(buffer);
+ cc = sc.toCartesian();
- sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
-
/* it is possible that if on phi=0 boundary, significant digits can push phi over 2pi */
- (E->trace.keep_within_bounds)(E,&x,&y,&z,&theta,&phi,&rad);
+ (E->trace.keep_within_bounds)(E,cc,sc);
Tracer new_tracer;
- new_tracer.setCoords(SphericalCoord(theta, phi, rad), CartesianCoord(x, y, z));
+ new_tracer.setCoords(sc, cc);
new_tracer.set_flavor(buffer[3]);
E->trace.tracers[j].push_back(new_tracer);
@@ -1125,9 +1093,6 @@
}
-
-
-
/*********** CHECK SUM **************************************************/
/* */
/* This functions checks to make sure number of tracers is preserved */
@@ -1192,7 +1157,6 @@
*theta=atan2(sqrt(temp),z);
*phi=myatan(y,x);
-
return;
}
@@ -1219,8 +1183,6 @@
return;
}
-
-
static void init_tracer_flavors(struct All_variables *E)
{
int j, kk, i;
@@ -1400,7 +1362,6 @@
int kk;
-
if (E->trace.nflavors > 0) {
E->trace.ntracer_flavor[j]=(int **)malloc(E->trace.nflavors*sizeof(int*));
for (kk=0;kk<E->trace.nflavors;kk++) {
@@ -1412,18 +1373,12 @@
}
}
-
fflush(E->trace.fpt);
return;
}
-
-
-
-
-
/********** ICHECK PROCESSOR SHELL *************/
/* returns -99 if rad is below current shell */
/* returns 0 if rad is above current shell */
@@ -1452,7 +1407,6 @@
if (rad<bottom_r) return -99;
-
/* Check top */
if (rad<top_r) return 1;
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h 2011-07-22 20:35:08 UTC (rev 18797)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/prototypes.h 2011-07-22 21:59:35 UTC (rev 18798)
@@ -152,10 +152,10 @@
void full_lost_souls(struct All_variables *);
void full_get_shape_functions(struct All_variables *, double [9], int, double, double, double);
double full_interpolate_data(struct All_variables *, double [9], double [9]);
-void full_get_velocity(struct All_variables *, int, int, double, double, double, double *);
-int full_icheck_cap(struct All_variables *, int, double, double, double, double);
-int full_iget_element(struct All_variables *, int, int, double, double, double, double, double, double);
-void full_keep_within_bounds(struct All_variables *, double *, double *, double *, double *, double *, double *);
+CartesianCoord full_get_velocity(struct All_variables *, int, int, SphericalCoord);
+int full_icheck_cap(struct All_variables *, int, CartesianCoord, double);
+int full_iget_element(struct All_variables *, int, int, CartesianCoord, SphericalCoord);
+void full_keep_within_bounds(struct All_variables *, CartesianCoord &, SphericalCoord &);
void analytical_test(struct All_variables *);
void analytical_runge_kutte(struct All_variables *, int, double, double *, double *, double *, double *, double *);
void analytical_test_function(struct All_variables *, double, double, double, double *, double *);
@@ -422,14 +422,14 @@
void regional_coord_of_cap(struct All_variables *, int, int);
/* Regional_tracer_advection.c */
void regional_tracer_setup(struct All_variables *);
-int regional_iget_element(struct All_variables *, int, int, double, double, double, double, double, double);
+int regional_iget_element(struct All_variables *, int, int, CartesianCoord, SphericalCoord);
int isearch_all(double *, int, double);
int isearch_neighbors(double *, int, double, int);
-int regional_icheck_cap(struct All_variables *, int, double, double, double, double);
+int regional_icheck_cap(struct All_variables *, int, SphericalCoord, double);
void regional_get_shape_functions(struct All_variables *, double [9], int, double, double, double);
double regional_interpolate_data(struct All_variables *, double [9], double [9]);
-void regional_get_velocity(struct All_variables *, int, int, double, double, double, double *);
-void regional_keep_within_bounds(struct All_variables *, double *, double *, double *, double *, double *, double *);
+CartesianCoord regional_get_velocity(struct All_variables *, int, int, SphericalCoord);
+void regional_keep_within_bounds(struct All_variables *, CartesianCoord &, SphericalCoord &);
void regional_lost_souls(struct All_variables *);
/* Regional_version_dependent.c */
void regional_node_locations(struct All_variables *);
Modified: mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h 2011-07-22 20:35:08 UTC (rev 18797)
+++ mc/3D/CitcomS/branches/eheien_dev/lib/tracer_defs.h 2011-07-22 21:59:35 UTC (rev 18798)
@@ -253,12 +253,11 @@
/*********************/
int (* iget_element)(struct All_variables*, int, int,
- double, double, double, double, double, double);
+ CartesianCoord, SphericalCoord);
- void (* get_velocity)(struct All_variables*, int, int,
- double, double, double, double*);
+ CartesianCoord (* get_velocity)(struct All_variables*, int, int,
+ SphericalCoord);
void (* keep_within_bounds)(struct All_variables*,
- double*, double*, double*,
- double*, double*, double*);
+ CartesianCoord &, SphericalCoord &);
};
More information about the CIG-COMMITS
mailing list