https://github.com/geodynamics/citcoms
Tip revision: b5721c162b1b01c9fdbfb544f8d60bf3742dd916 authored by Eric Heien on 02 August 2011, 16:52:33 UTC
Further work in converting tracer code to C++
Further work in converting tracer code to C++
Tip revision: b5721c1
Full_tracer_advection.c
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
* Clint Conrad, Michael Gurnis, and Eun-seo Choi.
* Copyright (C) 1994-2005, California Institute of Technology.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <math.h>
#include "element_definitions.h"
#include "global_defs.h"
#include "parsing.h"
#include "parallel_related.h"
#include "composition_related.h"
static void get_2dshape(struct All_variables *E,
int j, ElementID nelem,
CoordUV uv,
int iwedge, double * shape2d);
static void get_radial_shape(struct All_variables *E,
int j, ElementID nelem,
double rad, double *shaperad);
static CoordUV spherical_to_uv(struct All_variables *E,
SphericalCoord sc);
static int icheck_column_neighbors(struct All_variables *E,
int j, ElementID nel,
CartesianCoord cc,
double rad);
static int icheck_all_columns(struct All_variables *E,
int j,
CartesianCoord cc,
double rad);
static int icheck_element(struct All_variables *E,
int j, ElementID nel,
CartesianCoord cc,
double rad);
static int icheck_shell(struct All_variables *E,
ElementID nel, double rad);
static int icheck_element_column(struct All_variables *E,
int j, ElementID nel,
CartesianCoord cc,
double rad);
static int icheck_bounds(struct All_variables *E,
CartesianCoord test_point,
const CapBoundary bounds);
static double findradial(CartesianCoord vec,
BoundaryPoint bp);
static void fix_angle(double *angle);
static int iget_radial_element(struct All_variables *E,
int j, int iel,
double rad);
static int iget_regel(struct All_variables *E, int j,
double theta, double phi,
int *ntheta, int *nphi);
static void full_put_lost_tracers(struct All_variables *E,
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,
CartesianCoord cc, double rad);
/******* FULL TRACER INPUT *********************/
void full_tracer_input(struct All_variables *E)
{
int m = E->parallel.me;
/* Regular grid parameters */
/* (first fill uniform del[0] value) */
/* (later, in make_regular_grid, will adjust and distribute to caps */
E->trace.deltheta[0]=1.0;
E->trace.delphi[0]=1.0;
input_double("regular_grid_deltheta",&(E->trace.deltheta[0]),"1.0",m);
input_double("regular_grid_delphi",&(E->trace.delphi[0]),"1.0",m);
/* Analytical Test Function */
E->trace.ianalytical_tracer_test=0;
/* input_int("analytical_tracer_test",&(E->trace.ianalytical_tracer_test),
"0,0,nomax",m); */
}
/************** LOST SOULS ***************************************************/
/* */
/* This function is used to transport tracers to proper processor domains. */
/* (MPI parallel) */
/* All of the tracers that were sent to rlater arrays are destined to another*/
/* cap and sent there. Then they are raised up or down for multiple z procs. */
/* isend[j][n]=number of tracers this processor cap is sending to cap n */
/* ireceive[j][n]=number of tracers this processor cap receiving from cap n */
void full_lost_souls(struct All_variables *E)
{
/* This code works only if E->sphere.caps_per_proc==1 */
const int j = 1;
int ithiscap;
int ithatcap=1;
int isend[13][13];
int ireceive[13][13];
int isize[13];
int kk,pp;
int mm;
int numtracers;
int icheck=0;
int isend_position;
int ipos,ipos2,ipos3;
int idb;
int idestination_proc=0;
int isource_proc;
int isend_z[13][3];
int ireceive_z[13][3];
int isum[13];
int irad;
int ival;
int ithat_processor;
int ireceive_position;
int ivertical_neighbor;
int ilast_receiver_position;
int it;
int irec[13];
int irec_position;
int iel;
int num_tracers;
int isize_send;
int isize_receive;
int itemp_size;
int itracers_subject_to_vertical_transport[13];
double x,y,z;
double theta,phi,rad;
double *send[13][13];
double *receive[13][13];
double *send_z[13][3];
double *receive_z[13][3];
double *REC[13];
double begin_time = CPU_time0();
int number_of_caps=12;
int lev=E->mesh.levmax;
int num_ngb = E->parallel.TNUM_PASS[lev][j];
Tracer temp_tracer;
/* Note, if for some reason, the number of neighbors exceeds */
/* 50, which is unlikely, the MPI arrays must be increased. */
MPI_Status status[200];
MPI_Request request[200];
MPI_Status status1;
MPI_Status status2;
int itag=1;
parallel_process_sync(E);
if(E->control.verbose)
fprintf(E->trace.fpt, "Entering lost_souls()\n");
E->trace.istat_isend=E->trace.escaped_tracers[j].size();
/** debug **
for (kk=1; kk<=E->trace.istat_isend; kk++) {
fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
E->trace.rlater[j][0][kk],
E->trace.rlater[j][1][kk],
E->trace.rlater[j][2][kk]);
}
fflush(E->trace.fpt);
/**/
/* initialize isend and ireceive */
/* # of neighbors in the horizontal plane */
isize[j]=E->trace.escaped_tracers[j].size()*temp_tracer.size();
for (kk=0;kk<=num_ngb;kk++) isend[j][kk]=0;
for (kk=0;kk<=num_ngb;kk++) ireceive[j][kk]=0;
/* Allocate Maximum Memory to Send Arrays */
itemp_size=citmax(isize[j],1);
for (kk=0;kk<=num_ngb;kk++) {
if ((send[j][kk]=(double *)malloc(itemp_size*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
fflush(E->trace.fpt);
exit(10);
}
}
/** debug **
ithiscap=E->sphere.capid[j];
for (kk=1;kk<=num_ngb;kk++) {
ithatcap=E->parallel.PROCESSOR[lev][j].pass[kk];
fprintf(E->trace.fpt,"cap: %d me %d TNUM: %d rank: %d\n",
ithiscap,E->parallel.me,kk,ithatcap);
}
fflush(E->trace.fpt);
/**/
/* Pre communication */
full_put_lost_tracers(E, isend, send);
/* Send info to other processors regarding number of send tracers */
/* idb is the request array index variable */
/* Each send and receive has a request variable */
idb=0;
ithiscap=0;
/* if tracer is in same cap (nprocz>1) */
if (E->parallel.nprocz>1) {
ireceive[j][ithiscap]=isend[j][ithiscap];
}
for (kk=1;kk<=num_ngb;kk++) {
idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
MPI_Isend(&isend[j][kk],1,MPI_INT,idestination_proc,
11,E->parallel.world,&request[idb++]);
MPI_Irecv(&ireceive[j][kk],1,MPI_INT,idestination_proc,
11,E->parallel.world,&request[idb++]);
} /* end kk, number of neighbors */
/* Wait for non-blocking calls to complete */
MPI_Waitall(idb,request,status);
/** debug **
for (kk=0;kk<=num_ngb;kk++) {
if(kk==0)
isource_proc=E->parallel.me;
else
isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
fprintf(E->trace.fpt,"%d send %d to proc %d\n",
E->parallel.me,isend[j][kk],isource_proc);
fprintf(E->trace.fpt,"%d recv %d from proc %d\n",
E->parallel.me,ireceive[j][kk],isource_proc);
}
/**/
/* Allocate memory in receive arrays */
for (ithatcap=0;ithatcap<=num_ngb;ithatcap++) {
isize[j]=ireceive[j][ithatcap]*temp_tracer.size();
itemp_size=citmax(1,isize[j]);
if ((receive[j][ithatcap]=(double *)malloc(itemp_size*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
fflush(E->trace.fpt);
exit(10);
}
}
/* Now, send the tracers to proper caps */
idb=0;
ithiscap=0;
/* same cap */
if (E->parallel.nprocz>1) {
ithatcap=ithiscap;
isize[j]=isend[j][ithatcap]*temp_tracer.size();
for (mm=0;mm<isize[j];mm++) {
receive[j][ithatcap][mm]=send[j][ithatcap][mm];
}
}
/* neighbor caps */
for (kk=1;kk<=num_ngb;kk++) {
idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
isize[j]=isend[j][kk]*temp_tracer.size();
MPI_Isend(send[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
isize[j]=ireceive[j][kk]*temp_tracer.size();
MPI_Irecv(receive[j][kk],isize[j],MPI_DOUBLE,idestination_proc,
11,E->parallel.world,&request[idb++]);
} /* end kk, number of neighbors */
/* Wait for non-blocking calls to complete */
MPI_Waitall(idb,request,status);
/* Put all received tracers in array REC[j] */
/* This makes things more convenient. */
/* Sum up size of receive arrays (all tracers sent to this processor) */
isum[j]=0;
ithiscap=0;
for (kk=0;kk<=num_ngb;kk++) {
isum[j]=isum[j]+ireceive[j][kk];
}
itracers_subject_to_vertical_transport[j]=isum[j];
/* Allocate Memory for REC array */
isize[j]=isum[j]*temp_tracer.size();
isize[j]=citmax(isize[j],1);
if ((REC[j]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (g323)\n");
fflush(E->trace.fpt);
exit(10);
}
/* Put Received tracers in REC */
irec[j]=0;
irec_position=0;
for (kk=0;kk<=num_ngb;kk++) {
ithatcap=kk;
for (pp=0;pp<ireceive[j][ithatcap];pp++) {
irec[j]++;
ipos=pp*temp_tracer.size();
for (mm=0;mm<temp_tracer.size();mm++) {
ipos2=ipos+mm;
REC[j][irec_position]=receive[j][ithatcap][ipos2];
irec_position++;
} /* end mm (cycling tracer quantities) */
} /* end pp (cycling tracers) */
} /* end kk (cycling neighbors) */
/* Done filling REC */
/* VERTICAL COMMUNICATION */
if (E->parallel.nprocz>1) {
/* Allocate memory for send_z */
/* Make send_z the size of receive array (max size) */
/* (No dynamic reallocation of send_z necessary) */
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
isize[j]=itracers_subject_to_vertical_transport[j]*temp_tracer.size();
isize[j]=citmax(isize[j],1);
if ((send_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (c721)\n");
fflush(E->trace.fpt);
exit(10);
}
}
for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
ithat_processor=E->parallel.PROCESSORz[lev].pass[ivertical_neighbor];
/* initialize isend_z and ireceive_z array */
isend_z[j][ivertical_neighbor]=0;
ireceive_z[j][ivertical_neighbor]=0;
/* sort through receive array and check radius */
it=0;
num_tracers=irec[j];
for (kk=1;kk<=num_tracers;kk++) {
ireceive_position=it*temp_tracer.size();
it++;
irad=ireceive_position+2;
rad=REC[j][irad];
ival=icheck_that_processor_shell(E,j,ithat_processor,rad);
/* if tracer is in other shell, take out of receive array and give to send_z*/
if (ival==1) {
isend_position=isend_z[j][ivertical_neighbor]*temp_tracer.size();
isend_z[j][ivertical_neighbor]++;
ilast_receiver_position=(irec[j]-1)*temp_tracer.size();
for (mm=0;mm<temp_tracer.size();mm++) {
ipos=ireceive_position+mm;
ipos2=isend_position+mm;
send_z[j][ivertical_neighbor][ipos2]=REC[j][ipos];
/* eject tracer info from REC array, and replace with last tracer in array */
ipos3=ilast_receiver_position+mm;
REC[j][ipos]=REC[j][ipos3];
}
it--;
irec[j]--;
} /* end if ival===1 */
/* Otherwise, leave tracer */
} /* end kk (cycling through tracers) */
} /* end ivertical_neighbor */
/* Send arrays are now filled. */
/* Now send send information to vertical processor neighbor */
idb=0;
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
idestination_proc = E->parallel.PROCESSORz[lev].pass[kk];
MPI_Isend(&isend_z[j][kk],1,MPI_INT,idestination_proc,
14,E->parallel.world,&request[idb++]);
MPI_Irecv(&ireceive_z[j][kk],1,MPI_INT,idestination_proc,
14,E->parallel.world,&request[idb++]);
} /* end ivertical_neighbor */
/* Wait for non-blocking calls to complete */
MPI_Waitall(idb,request,status);
/** debug **
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
fprintf(E->trace.fpt, "PROC: %d IVN: %d (P: %d) "
"SEND: %d REC: %d\n",
E->parallel.me,kk,E->parallel.PROCESSORz[lev].pass[kk],
isend_z[j][kk],ireceive_z[j][kk]);
}
fflush(E->trace.fpt);
/**/
/* Allocate memory to receive_z arrays */
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
isize[j]=ireceive_z[j][kk]*temp_tracer.size();
isize[j]=citmax(isize[j],1);
if ((receive_z[j][kk]=(double *)malloc(isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (t590)\n");
fflush(E->trace.fpt);
exit(10);
}
}
/* Send Tracers */
idb=0;
for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++) {
idestination_proc = E->parallel.PROCESSORz[lev].pass[kk];
isize_send=isend_z[j][kk]*temp_tracer.size();
MPI_Isend(send_z[j][kk],isize_send,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
isize_receive=ireceive_z[j][kk]*temp_tracer.size();
MPI_Irecv(receive_z[j][kk],isize_receive,MPI_DOUBLE,idestination_proc,
15,E->parallel.world,&request[idb++]);
}
/* Wait for non-blocking calls to complete */
MPI_Waitall(idb,request,status);
/* Put tracers into REC array */
/* First, reallocate memory to REC */
isum[j]=0;
for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
isum[j] += ireceive_z[j][ivertical_neighbor];
}
isum[j] += irec[j];
isize[j]=isum[j]*temp_tracer.size();
if (isize[j]>0) {
if ((REC[j]=(double *)realloc(REC[j],isize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (i981)\n");
fprintf(E->trace.fpt,"isize: %d\n",isize[j]);
fflush(E->trace.fpt);
exit(10);
}
}
for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
for (kk=0;kk<ireceive_z[j][ivertical_neighbor];kk++) {
irec_position=irec[j]*temp_tracer.size();
irec[j]++;
ireceive_position=kk*temp_tracer.size();
for (mm=0;mm<temp_tracer.size();mm++) {
REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
}
}
}
/* Free Vertical Arrays */
for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++) {
free(send_z[j][ivertical_neighbor]);
free(receive_z[j][ivertical_neighbor]);
}
} /* endif nprocz>1 */
/* END OF VERTICAL TRANSPORT */
/* Put away tracers */
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]);
sc = new_tracer.getSphericalPos();
cc = new_tracer.getCartesianPos();
iel=(E->trace.iget_element)(E,j,UNDEFINED_ELEMENT,cc,sc);
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",
cc._x,cc._y,cc._z,sc._theta,sc._phi,sc._rad);
fflush(E->trace.fpt);
exit(10);
}
new_tracer.set_ielement(iel);
E->trace.tracers[j].push_back(new_tracer);
}
if(E->control.verbose){
fprintf(E->trace.fpt,"Freeing memory in lost_souls()\n");
fflush(E->trace.fpt);
}
parallel_process_sync(E);
/* Free Arrays */
free(REC[j]);
for (kk=0;kk<=num_ngb;kk++) {
free(send[j][kk]);
free(receive[j][kk]);
}
if(E->control.verbose){
fprintf(E->trace.fpt,"Leaving lost_souls()\n");
fflush(E->trace.fpt);
}
E->trace.lost_souls_time += CPU_time0() - begin_time;
}
static void full_put_lost_tracers(struct All_variables *E,
int isend[13][13], double *send[13][13])
{
const int j = 1;
int kk, pp;
int ithatcap, icheck;
int isend_position, ipos;
int lev = E->mesh.levmax;
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) {
cc = tr->getCartesianPos();
/* first check same cap if nprocz>1 */
if (E->parallel.nprocz>1) {
ithatcap=0;
icheck=full_icheck_cap(E,ithatcap,cc,tr->rad());
if (icheck==1) goto foundit;
}
/* check neighboring caps */
for (pp=1;pp<=E->parallel.TNUM_PASS[lev][j];pp++) {
ithatcap=pp;
icheck=full_icheck_cap(E,ithatcap,cc,tr->rad());
if (icheck==1) goto foundit;
}
/* 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",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);
exit(10);
}
foundit:
isend[j][ithatcap]++;
/* assign tracer to send */
isend_position=(isend[j][ithatcap]-1)*tr->size();
tr->writeToMem(&send[j][ithatcap][isend_position]);
} /* end kk, assigning tracers */
E->trace.escaped_tracers[j].clear();
}
/************************ GET SHAPE FUNCTION *********************************/
/* Real theta,phi,rad space is transformed into u,v space. This transformation */
/* maps great circles into straight lines. Here, elements boundaries are */
/* assumed to be great circle planes (not entirely true, it is actually only */
/* the nodal arrangement that lies upon great circles). Element boundaries */
/* are then mapped into planes. The element is then divided into 2 wedges */
/* in which standard shape functions are used to interpolate velocity. */
/* This transformation was found on the internet (refs were difficult to */
/* to obtain). It was tested that nodal configuration is indeed transformed */
/* into straight lines. */
/* Radial and azimuthal shape functions are decoupled. First find the shape */
/* functions associated with the 2D surface plane, then apply radial shape */
/* functions. */
/* */
/* Wedge information: */
/* */
/* Wedge 1 Wedge 2 */
/* _______ _______ */
/* */
/* wedge_node real_node wedge_node real_node */
/* ---------- --------- ---------- --------- */
/* */
/* 1 1 1 1 */
/* 2 2 2 3 */
/* 3 3 3 4 */
/* 4 5 4 5 */
/* 5 6 5 7 */
/* 6 7 6 8 */
void full_get_shape_functions(struct All_variables *E, int &shape_iwedge,
double shp[6], ElementID nelem,
SphericalCoord sc)
{
const int j = 1;
int iwedge,inum;
int i, kk;
int ival;
int itry;
CoordUV uv;
double shape2d[3];
double shaperad[2];
double shape[7];
int maxlevel=E->mesh.levmax;
const double eps=-1e-4;
/* find u and v using spherical coordinates */
uv = spherical_to_uv(E,sc);
inum=0;
itry=1;
try_again:
/* Check first wedge (1 of 2) */
iwedge=0;
next_wedge:
/* determine shape functions of wedge */
/* There are 3 shape functions for the triangular wedge */
get_2dshape(E,j,nelem,uv,iwedge,shape2d);
/* if any shape functions are negative, goto next wedge */
if (shape2d[0]<eps||shape2d[1]<eps||shape2d[2]<eps)
{
inum++;
/* AKMA clean this up */
if (inum>3)
{
fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>3!\n");
fflush(E->trace.fpt);
exit(10);
}
if (inum>1 && itry==1)
{
CartesianCoord cc;
fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>1\n");
fprintf(E->trace.fpt,"shape %f %f %f\n",shape2d[0],shape2d[1],shape2d[2]);
fprintf(E->trace.fpt,"u %f v %f element: %d \n",uv.u,uv.v, nelem);
fprintf(E->trace.fpt,"Element uv boundaries: \n");
for(kk=1;kk<=4;kk++) {
i = (E->ien[j][nelem].node[kk] - 1) / E->lmesh.noz + 1;
fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,(*E->trace.gnomonic)[i].u,(*E->trace.gnomonic)[i].v);
}
fprintf(E->trace.fpt,"theta: %f phi: %f rad: %f\n",sc._theta,sc._phi,sc._rad);
fprintf(E->trace.fpt,"Element theta-phi boundaries: \n");
for(kk=1;kk<=4;kk++)
fprintf(E->trace.fpt,"%d: Theta: %f Phi:%f\n",kk,E->sx[j][1][E->ien[j][nelem].node[kk]],E->sx[j][2][E->ien[j][nelem].node[kk]]);
cc = sc.toCartesian();
ival=icheck_element(E,j,nelem,cc,sc._rad);
fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
ival=(E->trace.iget_element)(E, j, UNDEFINED_ELEMENT, cc, sc);
fprintf(E->trace.fpt,"New Element?: %d\n",ival);
ival=icheck_column_neighbors(E,j,nelem,cc,sc._rad);
fprintf(E->trace.fpt,"New Element (neighs)?: %d\n",ival);
nelem=ival;
ival=icheck_element(E,j,nelem,cc,sc._rad);
fprintf(E->trace.fpt,"ICHECK?: %d\n",ival);
itry++;
if (ival>0) goto try_again;
fprintf(E->trace.fpt,"NO LUCK\n");
fflush(E->trace.fpt);
exit(10);
}
iwedge=1;
goto next_wedge;
}
/* Determine radial shape functions */
/* There are 2 shape functions radially */
get_radial_shape(E,j,nelem,sc._rad,shaperad);
/* There are 6 nodes to the solid wedge. */
/* The 6 shape functions assocated with the 6 nodes */
/* are products of radial and wedge shape functions. */
/* Sum of shape functions is 1 */
shape_iwedge = iwedge;
shp[0]=shaperad[0]*shape2d[0];
shp[1]=shaperad[0]*shape2d[1];
shp[2]=shaperad[0]*shape2d[2];
shp[3]=shaperad[1]*shape2d[0];
shp[4]=shaperad[1]*shape2d[1];
shp[5]=shaperad[1]*shape2d[2];
/** debug **
fprintf(E->trace.fpt, "shp: %e %e %e %e %e %e\n",
shp[0], shp[1], shp[2], shp[3], shp[4], shp[5]);
/**/
}
/************************ GET VELOCITY ***************************************/
/* */
/* This function interpolates tracer velocity using gnominic interpolation. */
/* The element is divided into 2 wedges in which standard shape functions */
/* are used to interpolate velocity. */
/* */
/* Wedge information: */
/* */
/* Wedge 1 Wedge 2 */
/* _______ _______ */
/* */
/* wedge_node real_node wedge_node real_node */
/* ---------- --------- ---------- --------- */
/* */
/* 1 1 1 1 */
/* 2 2 2 3 */
/* 3 3 3 4 */
/* 4 5 4 5 */
/* 5 6 5 7 */
/* 6 7 6 8 */
CartesianCoord full_get_velocity(struct All_variables *E,
int j, ElementID nelem,
SphericalCoord sc)
{
int iwedge, i;
const int sphere_key = 0;
double shape[6];
CartesianCoord VV[9], vel[6], result;
full_get_shape_functions(E, iwedge, shape, nelem, sc);
/* get cartesian velocity */
velo_from_element_d(E, VV, j, nelem, sphere_key);
/* depending on wedge, set up velocity points */
if (iwedge==0) {
vel[0]=VV[1];
vel[1]=VV[2];
vel[2]=VV[3];
vel[3]=VV[5];
vel[4]=VV[6];
vel[5]=VV[7];
} else if (iwedge==1) {
vel[0]=VV[1];
vel[1]=VV[3];
vel[2]=VV[4];
vel[3]=VV[5];
vel[4]=VV[7];
vel[5]=VV[8];
}
for (i=0;i<6;++i) result += vel[i] * shape[i];
return result;
}
/***************************************************************/
/* GET 2DSHAPE */
/* */
/* This function determines shape functions at u,v */
/* This method uses standard linear shape functions of */
/* triangular elements. (See Cuvelier, Segal, and */
/* van Steenhoven, 1986). */
static void get_2dshape(struct All_variables *E,
int j, ElementID nelem,
CoordUV uv,
int iwedge, double * shape2d)
{
/* convert nelem to surface element number */
int n = (nelem - 1) / E->lmesh.elz + 1;
/* shape functions */
shape2d[0] = E->trace.shape_coefs[j][iwedge][n]->applyShapeFunc(0, uv);
shape2d[1] = E->trace.shape_coefs[j][iwedge][n]->applyShapeFunc(1, uv);
shape2d[2] = E->trace.shape_coefs[j][iwedge][n]->applyShapeFunc(2, uv);
/** debug **
fprintf(E->trace.fpt, "el=%d els=%d iwedge=%d shape=(%e %e %e)\n",
nelem, n, iwedge, shape2d[0], shape2d[1], shape2d[2]);
/**/
}
/***************************************************************/
/* GET RADIAL SHAPE */
/* */
/* This function determines radial shape functions at rad */
static void get_radial_shape(struct All_variables *E,
int j, ElementID nelem,
double rad, double *shaperad)
{
int node1,node5;
double rad1,rad5,f1,f2,delrad;
const double eps=1e-6;
double top_bound=1.0+eps;
double bottom_bound=0.0-eps;
node1=E->ien[j][nelem].node[1];
node5=E->ien[j][nelem].node[5];
rad1=E->sx[j][3][node1];
rad5=E->sx[j][3][node5];
delrad=rad5-rad1;
f1=(rad-rad1)/delrad;
f2=(rad5-rad)/delrad;
/* Save a small amount of computation here */
/* because f1+f2=1, shapes can be switched */
/*
shaperad[0]=1.0-f1=1.0-(1.0-f2)=f2;
shaperad[1]=1.0-f2=1.0-(10-f1)=f1;
*/
shaperad[0]=f2;
shaperad[1]=f1;
/* Some error control */
if (shaperad[0]>top_bound || shaperad[0]<bottom_bound ||
shaperad[1]>top_bound || shaperad[1]<bottom_bound)
{
fprintf(E->trace.fpt,"ERROR(get_radial_shape)\n");
fprintf(E->trace.fpt,"shaperad[0]: %f \n",shaperad[0]);
fprintf(E->trace.fpt,"shaperad[1]: %f \n",shaperad[1]);
fflush(E->trace.fpt);
exit(10);
}
}
/**************************************************************/
/* SPHERICAL TO UV */
/* */
/* This function transforms theta and phi to new coords */
/* u and v using gnomonic projection. */
static CoordUV spherical_to_uv(struct All_variables *E,
SphericalCoord sc)
{
double phi_f;
double cosc;
double cos_theta_f,sin_theta_f;
double cost,sint,cosp2,sinp2;
/* theta_f and phi_f are the reference points of the cap */
phi_f = E->trace.gnomonic_reference_phi;
cos_theta_f = (*E->trace.gnomonic)[0].u;
sin_theta_f = (*E->trace.gnomonic)[0].v;
cost=cos(sc._theta);
/*
sint=sin(theta);
*/
sint=sqrt(1.0-cost*cost);
cosp2=cos(sc._phi-phi_f);
sinp2=sin(sc._phi-phi_f);
cosc=cos_theta_f*cost+sin_theta_f*sint*cosp2;
cosc=1.0/cosc;
return CoordUV(sint*sinp2*cosc, (sin_theta_f*cost-cos_theta_f*sint*cosp2)*cosc);
/** debug **
fprintf(E->trace.fpt, "(%e %e) -> (%e %e)\n",
theta, phi, *u, *v);
/**/
}
/*********** MAKE REGULAR GRID ********************************/
/* */
/* This function generates the finer regular grid which is */
/* mapped to real elements */
void make_regular_grid(struct All_variables *E)
{
int j;
int kk;
int mm;
int pp,node;
int numtheta,numphi;
int nodestheta,nodesphi;
unsigned int numregel;
unsigned int numregnodes;
int idum1,idum2;
int ifound_one;
int ival;
int ilast_el;
int imap;
int elz;
int nelsurf;
int iregnode[5];
int ntheta,nphi;
int ichoice;
int icount;
int itemp[5];
int iregel;
int istat_ichoice[13][5];
int isum;
double x,y,z;
double theta,phi,rad;
double deltheta;
double delphi;
double thetamax,thetamin;
double phimax,phimin;
double start_time;
double theta_min,phi_min;
double theta_max,phi_max;
double half_diff;
double expansion;
double *tmin;
double *tmax;
double *fmin;
double *fmax;
const double two_pi=2.0*M_PI;
elz=E->lmesh.elz;
nelsurf=E->lmesh.elx*E->lmesh.ely;
//TODO: find the bounding box of the mesh, if the box is too close to
// to core, set a flag (rotated_reggrid) to true and rotate the
// bounding box to the equator. Generate the regular grid with the new
// bounding box. The rotation should be a simple one, e.g.
// (theta, phi) -> (??)
// Whenever the regular grid is used, check the flat (rotated_reggrid),
// if true, rotate the checkpoint as well.
/* note, mesh is rotated along theta 22.5 degrees divided by elx. */
/* We at least want that much expansion here! Otherwise, theta min */
/* will not be valid near poles. We do a little more (x2) to be safe */
/* Typically 1-2 degrees. Look in nodal_mesh.c for this. */
expansion=2.0*0.5*(M_PI/4.0)/(1.0*E->lmesh.elx);
start_time=CPU_time0();
if (E->parallel.me==0) fprintf(stderr,"Generating Regular Grid\n");
/* for each cap, determine theta and phi bounds, watch out near poles */
numregnodes=0;
for(j=1;j<=E->sphere.caps_per_proc;j++)
{
thetamax=0.0;
thetamin=M_PI;
phimax=two_pi;
phimin=0.0;
for (kk=1;kk<=E->lmesh.nno;kk=kk+E->lmesh.noz)
{
theta=E->sx[j][1][kk];
phi=E->sx[j][2][kk];
thetamax=citmax(thetamax,theta);
thetamin=citmin(thetamin,theta);
}
/* expand range slightly (should take care of poles) */
thetamax=thetamax+expansion;
thetamax=citmin(thetamax,M_PI);
thetamin=thetamin-expansion;
thetamin=citmax(thetamin,0.0);
/* Convert input data from degrees to radians */
deltheta=E->trace.deltheta[0]*M_PI/180.0;
delphi=E->trace.delphi[0]*M_PI/180.0;
/* Adjust deltheta and delphi to fit a uniform number of regular elements */
numtheta=fabs(thetamax-thetamin)/deltheta;
numphi=fabs(phimax-phimin)/delphi;
nodestheta=numtheta+1;
nodesphi=numphi+1;
numregel=numtheta*numphi;
numregnodes=nodestheta*nodesphi;
if ((numtheta==0)||(numphi==0))
{
fprintf(E->trace.fpt,"Error(make_regular_grid): numtheta: %d numphi: %d\n",numtheta,numphi);
fflush(E->trace.fpt);
exit(10);
}
deltheta=fabs(thetamax-thetamin)/(1.0*numtheta);
delphi=fabs(phimax-phimin)/(1.0*numphi);
/* fill global variables */
E->trace.deltheta[j]=deltheta;
E->trace.delphi[j]=delphi;
E->trace.numtheta[j]=numtheta;
E->trace.numphi[j]=numphi;
E->trace.thetamax[j]=thetamax;
E->trace.thetamin[j]=thetamin;
E->trace.phimax[j]=phimax;
E->trace.phimin[j]=phimin;
E->trace.numregel[j]=numregel;
E->trace.numregnodes[j]=numregnodes;
if ( ((1.0*numregel)/(1.0*E->lmesh.elx*E->lmesh.ely)) < 0.5 )
{
fprintf(E->trace.fpt,"\n ! WARNING: regular/real ratio low: %f ! \n",
((1.0*numregel)/(1.0*E->lmesh.nel)) );
fprintf(E->trace.fpt," Should reduce size of regular mesh\n");
fprintf(stderr,"! WARNING: regular/real ratio low: %f ! \n",
((1.0*numregel)/(1.0*E->lmesh.nel)) );
fprintf(stderr," Should reduce size of regular mesh\n");
fflush(E->trace.fpt);
if (E->trace.itracer_warnings) exit(10);
}
/* print some output */
fprintf(E->trace.fpt,"\nRegular grid:\n");
fprintf(E->trace.fpt,"Theta min: %f max: %f \n",thetamin,thetamax);
fprintf(E->trace.fpt,"Phi min: %f max: %f \n",phimin,phimax);
fprintf(E->trace.fpt,"Adjusted deltheta: %f delphi: %f\n",deltheta,delphi);
fprintf(E->trace.fpt,"(numtheta: %d numphi: %d)\n",numtheta,numphi);
fprintf(E->trace.fpt,"Number of regular elements: %d (nodes: %d)\n",numregel,numregnodes);
fprintf(E->trace.fpt,"regular/real ratio: %f\n",((1.0*numregel)/(1.0*E->lmesh.elx*E->lmesh.ely)));
fflush(E->trace.fpt);
/* Allocate memory for regnodetoel */
/* Regtoel is an integer array which represents nodes on */
/* the regular mesh. Each node on the regular mesh contains */
/* the real element value if one exists (UNDEFINED_ELEMENT otherwise) */
if ((E->trace.regnodetoel[j]=(int *)malloc((numregnodes+1)*sizeof(int)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular) -no memory - uh3ud\n");
fflush(E->trace.fpt);
exit(10);
}
/* Initialize regnodetoel - reg elements not used = UNDEFINED_ELEMENT */
for (kk=1;kk<=numregnodes;kk++)
{
E->trace.regnodetoel[j][kk]=UNDEFINED_ELEMENT;
}
/* Begin Mapping (only need to use surface elements) */
parallel_process_sync(E);
if (E->parallel.me==0) fprintf(stderr,"Beginning Mapping\n");
/* Generate temporary arrays of max and min values for each surface element */
if ((tmin=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
fflush(E->trace.fpt);
exit(10);
}
if ((tmax=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
fflush(E->trace.fpt);
exit(10);
}
if ((fmin=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
fflush(E->trace.fpt);
exit(10);
}
if ((fmax=(double *)malloc((nelsurf+1)*sizeof(double)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular) -no memory - 7t1a\n");
fflush(E->trace.fpt);
exit(10);
}
for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
{
kk=mm/elz;
theta_min=M_PI;
theta_max=0.0;
phi_min=two_pi;
phi_max=0.0;
for (pp=1;pp<=4;pp++)
{
node=E->ien[j][mm].node[pp];
theta=E->sx[j][1][node];
phi=E->sx[j][2][node];
theta_min=citmin(theta_min,theta);
theta_max=citmax(theta_max,theta);
phi_min=citmin(phi_min,phi);
phi_max=citmax(phi_max,phi);
}
/* add half difference to phi and expansion to theta to be safe */
theta_max=theta_max+expansion;
theta_min=theta_min-expansion;
theta_max=citmin(M_PI,theta_max);
theta_min=citmax(0.0,theta_min);
half_diff=0.5*(phi_max-phi_min);
phi_max=phi_max+half_diff;
phi_min=phi_min-half_diff;
fix_angle(&phi_max);
fix_angle(&phi_min);
if (phi_min>phi_max)
{
phi_min=0.0;
phi_max=two_pi;
}
tmin[kk]=theta_min;
tmax[kk]=theta_max;
fmin[kk]=phi_min;
fmax[kk]=phi_max;
}
/* end looking through elements */
ifound_one=0;
rad=E->sphere.ro;
imap=0;
for (kk=1;kk<=numregnodes;kk++)
{
E->trace.regnodetoel[j][kk]=UNDEFINED_ELEMENT;
/* find theta and phi for a given regular node */
idum1=(kk-1)/(numtheta+1);
idum2=kk-1-idum1*(numtheta+1);
theta=thetamin+(1.0*idum2*deltheta);
phi=phimin+(1.0*idum1*delphi);
SphericalCoord sc(theta,phi,rad);
CartesianCoord cc;
cc = sc.toCartesian();
ilast_el=1;
/* if previous element not found yet, check all surface elements */
/*
if (ifound_one==0)
{
for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
{
pp=mm/elz;
if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
{
ival=icheck_element_column(E,j,mm,x,y,z,rad);
if (ival>0)
{
ilast_el=mm;
ifound_one++;
E->trace.regnodetoel[j][kk]=mm;
goto foundit;
}
}
}
goto foundit;
}
*/
/* first check previous element */
ival=icheck_element_column(E,j,ilast_el,cc,rad);
if (ival>0)
{
E->trace.regnodetoel[j][kk]=ilast_el;
goto foundit;
}
/* check neighbors */
ival=icheck_column_neighbors(E,j,ilast_el,cc,rad);
if (ival>0)
{
E->trace.regnodetoel[j][kk]=ival;
ilast_el=ival;
goto foundit;
}
/* check all */
for (mm=elz;mm<=E->lmesh.nel;mm=mm+elz)
{
pp=mm/elz;
if ( (theta>=tmin[pp]) && (theta<=tmax[pp]) && (phi>=fmin[pp]) && (phi<=fmax[pp]) )
{
ival=icheck_element_column(E,j,mm,cc,rad);
if (ival>0)
{
ilast_el=mm;
E->trace.regnodetoel[j][kk]=mm;
goto foundit;
}
}
}
foundit:
if (E->trace.regnodetoel[j][kk]>0) imap++;
} /* end all regular nodes (kk) */
fprintf(E->trace.fpt,"percentage mapped: %f\n", (1.0*imap)/(1.0*numregnodes)*100.0);
fflush(E->trace.fpt);
/* free temporary arrays */
free(tmin);
free(tmax);
free(fmin);
free(fmax);
} /* end j */
/* some error control */
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
for (kk=1;kk<=numregnodes;kk++)
{
if (E->trace.regnodetoel[j][kk]!=UNDEFINED_ELEMENT)
{
if ( (E->trace.regnodetoel[j][kk]<1)||(E->trace.regnodetoel[j][kk]>E->lmesh.nel) )
{
fprintf(stderr,"Error(make_regular_grid)-invalid element: %d\n",E->trace.regnodetoel[j][kk]);
fprintf(E->trace.fpt,"Error(make_regular_grid)-invalid element: %d\n",E->trace.regnodetoel[j][kk]);
fflush(E->trace.fpt);
fflush(stderr);
exit(10);
}
}
}
}
/* Now put regnodetoel information into regtoel */
if (E->parallel.me==0) fprintf(stderr,"Beginning Regtoel submapping \n");
/* AKMA decided it would be more efficient to have reg element choice array */
/* rather than reg node array as used before */
for(j=1;j<=E->sphere.caps_per_proc;j++)
{
/* initialize statistical counter */
for (pp=0;pp<=4;pp++) istat_ichoice[j][pp]=0;
/* Allocate memory for regtoel */
/* Regtoel consists of 4 positions for each regular element */
/* Position[0] lists the number of element choices (later */
/* referred to as ichoice), followed */
/* by the possible element choices. */
/* ex.) A regular element has 4 nodes. Each node resides in */
/* a real element. The number of real elements a regular */
/* element touches (one of its nodes are in) is ichoice. */
/* Special ichoice notes: */
/* ichoice=-1 all regular element nodes = UNDEFINED_ELEMENT (no elements) */
/* ichoice=0 all 4 corners within one element */
/* ichoice=1 one element choice (diff from ichoice 0 in */
/* that perhaps one reg node is in an element */
/* and the rest are not (UNDEFINED_ELEMENT). */
/* ichoice>1 Multiple elements to check */
numregel= E->trace.numregel[j];
for (pp=0;pp<=4;pp++)
{
if ((E->trace.regtoel[j][pp]=(int *)malloc((numregel+1)*sizeof(int)))==NULL)
{
fprintf(E->trace.fpt,"ERROR(make regular)-no memory 98d (%d %d %d)\n",pp,numregel,j);
fflush(E->trace.fpt);
exit(10);
}
}
numtheta=E->trace.numtheta[j];
numphi=E->trace.numphi[j];
for (nphi=1;nphi<=numphi;nphi++)
{
for (ntheta=1;ntheta<=numtheta;ntheta++)
{
iregel=ntheta+(nphi-1)*numtheta;
/* initialize regtoel (not necessary really) */
for (pp=0;pp<=4;pp++) E->trace.regtoel[j][pp][iregel]=-33;
if ( (iregel>numregel)||(iregel<1) )
{
fprintf(E->trace.fpt,"ERROR(make_regular_grid)-weird iregel: %d (max: %d)\n",iregel,numregel);
fflush(E->trace.fpt);
exit(10);
}
iregnode[1]=iregel+(nphi-1);
iregnode[2]=iregel+nphi;
iregnode[3]=iregel+nphi+E->trace.numtheta[j]+1;
iregnode[4]=iregel+nphi+E->trace.numtheta[j];
for (kk=1;kk<=4;kk++)
{
if ((iregnode[kk]<1)||(iregnode[kk]>numregnodes))
{
fprintf(E->trace.fpt,"ERROR(make regular)-bad regnode %d\n",iregnode[kk]);
fflush(E->trace.fpt);
exit(10);
}
if (E->trace.regnodetoel[j][iregnode[kk]]>E->lmesh.nel)
{
fprintf(E->trace.fpt,"AABB HERE %d %d %d %d\n",iregel,iregnode[kk],kk,E->trace.regnodetoel[j][iregnode[kk]]);
fflush(E->trace.fpt);
}
}
/* find number of choices */
ichoice=0;
icount=0;
for (kk=1;kk<=4;kk++)
{
if (E->trace.regnodetoel[j][iregnode[kk]]<=0) goto next_corner;
icount++;
for (pp=1;pp<=(kk-1);pp++)
{
if (E->trace.regnodetoel[j][iregnode[kk]]==E->trace.regnodetoel[j][iregnode[pp]]) goto next_corner;
}
ichoice++;
itemp[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
if ((ichoice<0) || (ichoice>4) )
{
fprintf(E->trace.fpt,"ERROR(make regular) - weird ichoice %d \n",ichoice);
fflush(E->trace.fpt);
exit(10);
}
if ((itemp[ichoice]<0) || (itemp[ichoice]>E->lmesh.nel) )
{
fprintf(E->trace.fpt,"ERROR(make regular) - weird element choice %d %d\n",itemp[ichoice],ichoice);
fflush(E->trace.fpt);
exit(10);
}
next_corner:
;
} /* end kk */
istat_ichoice[j][ichoice]++;
if ((ichoice<0) || (ichoice>4))
{
fprintf(E->trace.fpt,"ERROR(make_regular)-wierd ichoice %d\n",ichoice);
fflush(E->trace.fpt);
exit(10);
}
if (ichoice==0)
{
E->trace.regtoel[j][0][iregel]=-1;
/*
fprintf(E->trace.fpt,"HH1: (%p) iregel: %d ichoice: %d value: %d %d\n",&E->trace.regtoel[j][1][iregel],iregel,ichoice,E->trace.regtoel[j][0][iregel],E->trace.regtoel[j][1][iregel]);
*/
}
else if ( (ichoice==1) && (icount==4) )
{
E->trace.regtoel[j][0][iregel]=0;
E->trace.regtoel[j][1][iregel]=itemp[1];
/*
fprintf(E->trace.fpt,"HH2: (%p) iregel: %d ichoice: %d value: %d %d\n",&E->trace.regtoel[j][1][iregel],iregel,ichoice,E->trace.regtoel[j][0][iregel],E->trace.regtoel[j][1][iregel]);
*/
if (itemp[1]<1 || itemp[1]>E->lmesh.nel)
{
fprintf(E->trace.fpt,"ERROR(make_regular)-huh? wierd itemp\n");
fflush(E->trace.fpt);
exit(10);
}
}
else if ( (ichoice>0) && (ichoice<5) )
{
E->trace.regtoel[j][0][iregel]=ichoice;
for (pp=1;pp<=ichoice;pp++)
{
E->trace.regtoel[j][pp][iregel]=itemp[pp];
/*
fprintf(E->trace.fpt,"HH:(%p) iregel: %d ichoice: %d pp: %d value: %d %d\n",&E->trace.regtoel[j][pp][iregel],iregel,ichoice,pp,itemp[pp],E->trace.regtoel[j][pp][iregel]);
*/
if (itemp[pp]<1 || itemp[pp]>E->lmesh.nel)
{
fprintf(E->trace.fpt,"ERROR(make_regular)-huh? wierd itemp 2 \n");
fflush(E->trace.fpt);
exit(10);
}
}
}
else
{
fprintf(E->trace.fpt,"ERROR(make_regular)- should not be here! %d\n",ichoice);
fflush(E->trace.fpt);
exit(10);
}
}
}
/* can now free regnodetoel */
free (E->trace.regnodetoel[j]);
/* testing */
for (kk=1;kk<=E->trace.numregel[j];kk++)
{
if ((E->trace.regtoel[j][0][kk]<-1)||(E->trace.regtoel[j][0][kk]>4))
{
fprintf(E->trace.fpt,"ERROR(make regular) regtoel ichoice0? %d %d \n",kk,E->trace.regtoel[j][pp][kk]);
fflush(E->trace.fpt);
exit(10);
}
for (pp=1;pp<=4;pp++)
{
if (((E->trace.regtoel[j][pp][kk]<1)&&(E->trace.regtoel[j][pp][kk]!=-33))||(E->trace.regtoel[j][pp][kk]>E->lmesh.nel))
{
fprintf(E->trace.fpt,"ERROR(make regular) (%p) regtoel? %d %d(%d) %d\n",&E->trace.regtoel[j][pp][kk],kk,pp,E->trace.regtoel[j][0][kk],E->trace.regtoel[j][pp][kk]);
fflush(E->trace.fpt);
exit(10);
}
}
}
} /* end j */
fprintf(E->trace.fpt,"Mapping completed (%f seconds)\n",CPU_time0()-start_time);
fflush(E->trace.fpt);
parallel_process_sync(E);
if (E->parallel.me==0) fprintf(stderr,"Mapping completed (%f seconds)\n",CPU_time0()-start_time);
/* Print out information regarding regular/real element coverage */
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
isum=0;
for (kk=0;kk<=4;kk++) isum=isum+istat_ichoice[j][kk];
fprintf(E->trace.fpt,"\n\nInformation regarding number of real elements per regular elements\n");
fprintf(E->trace.fpt," (stats done on regular elements that were used)\n");
fprintf(E->trace.fpt,"Ichoice is number of real elements touched by a regular element\n");
fprintf(E->trace.fpt," (ichoice=0 is optimal)\n");
fprintf(E->trace.fpt,"Ichoice=0: %f percent\n",(100.0*istat_ichoice[j][0])/(1.0*isum));
fprintf(E->trace.fpt,"Ichoice=1: %f percent\n",(100.0*istat_ichoice[j][1])/(1.0*isum));
fprintf(E->trace.fpt,"Ichoice=2: %f percent\n",(100.0*istat_ichoice[j][2])/(1.0*isum));
fprintf(E->trace.fpt,"Ichoice=3: %f percent\n",(100.0*istat_ichoice[j][3])/(1.0*isum));
fprintf(E->trace.fpt,"Ichoice=4: %f percent\n",(100.0*istat_ichoice[j][4])/(1.0*isum));
} /* end j */
}
/********* ICHECK COLUMN NEIGHBORS ***************************/
/* */
/* This function check whether a point is in a neighboring */
/* column. Neighbor surface element number is returned */
static int icheck_column_neighbors(struct All_variables *E,
int j, ElementID nel,
CartesianCoord cc,
double rad)
{
int ival;
int neighbor[25];
int elx,ely,elz;
int elxz;
int kk;
/*
const int number_of_neighbors=24;
*/
/* maybe faster to only check inner ring */
const int number_of_neighbors=8;
elx=E->lmesh.elx;
ely=E->lmesh.ely;
elz=E->lmesh.elz;
elxz=elx*elz;
/* inner ring */
neighbor[1]=nel-elxz-elz;
neighbor[2]=nel-elxz;
neighbor[3]=nel-elxz+elz;
neighbor[4]=nel-elz;
neighbor[5]=nel+elz;
neighbor[6]=nel+elxz-elz;
neighbor[7]=nel+elxz;
neighbor[8]=nel+elxz+elz;
/* outer ring */
neighbor[9]=nel+2*elxz-elz;
neighbor[10]=nel+2*elxz;
neighbor[11]=nel+2*elxz+elz;
neighbor[12]=nel+2*elxz+2*elz;
neighbor[13]=nel+elxz+2*elz;
neighbor[14]=nel+2*elz;
neighbor[15]=nel-elxz+2*elz;
neighbor[16]=nel-2*elxz+2*elz;
neighbor[17]=nel-2*elxz+elz;
neighbor[18]=nel-2*elxz;
neighbor[19]=nel-2*elxz-elz;
neighbor[20]=nel-2*elxz-2*elz;
neighbor[21]=nel-elxz-2*elz;
neighbor[22]=nel-2*elz;
neighbor[23]=nel+elxz-2*elz;
neighbor[24]=nel+2*elxz-2*elz;
for (kk=1;kk<=number_of_neighbors;kk++)
{
if ((neighbor[kk]>=1)&&(neighbor[kk]<=E->lmesh.nel))
{
ival=icheck_element_column(E,j,neighbor[kk],cc,rad);
if (ival>0)
{
return neighbor[kk];
}
}
}
return UNDEFINED_ELEMENT;
}
/********** ICHECK ALL COLUMNS ********************************/
/* */
/* This function check all columns until the proper one for */
/* a point (x,y,z) is found. The surface element is returned */
/* else UNDEFINED_ELEMENT if can't be found. */
static int icheck_all_columns(struct All_variables *E,
int j,
CartesianCoord cc,
double rad)
{
int icheck;
int nel;
int elz=E->lmesh.elz;
int numel=E->lmesh.nel;
for (nel=elz;nel<=numel;nel+=elz)
{
icheck=icheck_element_column(E,j,nel,cc,rad);
if (icheck==1) return nel;
}
return UNDEFINED_ELEMENT;
}
/******* ICHECK ELEMENT *************************************/
/* */
/* This function serves to determine if a point lies within */
/* a given element */
static int icheck_element(struct All_variables *E,
int j, ElementID nel,
CartesianCoord cc,
double rad)
{
int icheck;
icheck = icheck_shell(E, nel, rad);
if (icheck == 0) return 0;
icheck = icheck_element_column(E, j, nel, cc, rad);
if (icheck == 0) return 0;
return 1;
}
/******** ICHECK SHELL ************************************/
/* */
/* This function serves to check whether a point lies */
/* within the proper radial shell of a given element */
/* note: j set to 1; shouldn't depend on cap */
static int icheck_shell(struct All_variables *E,
ElementID nel, double rad)
{
int ibottom_node, itop_node;
double bottom_rad, top_rad;
ibottom_node=E->ien[1][nel].node[1];
itop_node=E->ien[1][nel].node[5];
bottom_rad=E->sx[1][3][ibottom_node];
top_rad=E->sx[1][3][itop_node];
if ((rad>=bottom_rad)&&(rad<top_rad)) return 1;
else return 0;
}
/******** ICHECK ELEMENT COLUMN ****************************/
/* */
/* This function serves to determine if a point lies within */
/* a given element's column */
static int icheck_element_column(struct All_variables *E,
int j, ElementID nel,
CartesianCoord cc,
double rad)
{
CartesianCoord test_point;
CapBoundary bounds;
int lev, kk, node;
lev = E->mesh.levmax;
E->trace.istat_elements_checked++;
/* surface coords of element nodes */
for (kk=0;kk<4;kk++)
{
node=E->ien[j][nel].node[kk+4+1];
bounds.setCartTrigBounds(kk,
CartesianCoord(E->x[j][1][node], E->x[j][2][node], E->x[j][3][node]),
E->SinCos[lev][j][2][node], /* cos(theta) */
E->SinCos[lev][j][0][node], /* sin(theta) */
E->SinCos[lev][j][3][node], /* cos(phi) */
E->SinCos[lev][j][1][node]); /* sin(phi) */
}
/* test_point - project to outer radius */
test_point = cc/rad;
return icheck_bounds(E,test_point,bounds);
}
/********* ICHECK CAP ***************************************/
/* */
/* This function serves to determine if a point lies within */
/* a given cap */
/* */
int full_icheck_cap(struct All_variables *E, int icap,
CartesianCoord cc, double rad)
{
CartesianCoord test_point;
/* test_point - project to outer radius */
test_point = cc/rad;
return icheck_bounds(E,test_point,E->trace.boundaries[icap]);
}
/***** ICHECK BOUNDS ******************************/
/* */
/* This function check if a test_point is bounded */
/* by 4 nodes */
/* This is done by: */
/* 1) generate vectors from node to node */
/* 2) generate vectors from each node to point */
/* in question */
/* 3) for each node, take cross product of vector */
/* pointing to it from previous node and */
/* vector from node to point in question */
/* 4) Find radial components of all the cross */
/* products. */
/* 5) If all radial components are positive, */
/* point is bounded by the 4 nodes */
/* 6) If some radial components are negative */
/* point is on a boundary - adjust it an */
/* epsilon amount for this analysis only */
/* which will force it to lie in one element */
/* or cap */
static int icheck_bounds(struct All_variables *E,
CartesianCoord test_point,
const CapBoundary bounds)
{
int number_of_tries=0;
CartesianCoord v12, v23, v34, v41;
CartesianCoord v1p, v2p, v3p, v4p;
CartesianCoord cross1, cross2, cross3, cross4;
SphericalCoord sc;
double rad1,rad2,rad3,rad4;
double tiny, eps;
/* make vectors from node to node */
v12 = bounds[1].cartesian() - bounds[0].cartesian();
v23 = bounds[2].cartesian() - bounds[1].cartesian();
v34 = bounds[3].cartesian() - bounds[2].cartesian();
v41 = bounds[0].cartesian() - bounds[3].cartesian();
try_again:
number_of_tries++;
/* make vectors from test point to node */
v1p = test_point - bounds[0].cartesian();
v2p = test_point - bounds[1].cartesian();
v3p = test_point - bounds[2].cartesian();
v4p = test_point - bounds[3].cartesian();
/* Calculate cross products */
cross2 = v12.crossProduct(v2p);
cross3 = v23.crossProduct(v3p);
cross4 = v34.crossProduct(v4p);
cross1 = v41.crossProduct(v1p);
/* Calculate radial component of cross products */
rad1=findradial(cross1,bounds[0]);
rad2=findradial(cross2,bounds[1]);
rad3=findradial(cross3,bounds[2]);
rad4=findradial(cross4,bounds[3]);
/* Check if any radial components is zero (along a boundary), adjust if so */
/* Hopefully, this doesn't happen often, may be expensive */
tiny=1e-15;
eps=1e-6;
if (number_of_tries>3)
{
fprintf(E->trace.fpt,"Error(icheck_bounds)-too many tries\n");
fprintf(E->trace.fpt,"Rads: %f %f %f %f\n",rad1,rad2,rad3,rad4);
fprintf(E->trace.fpt,"Test Point: %f %f %f \n",test_point._x,test_point._y,test_point._z);
//fprintf(E->trace.fpt,"Nodal points: 1: %f %f %f\n",rnode1._x,rnode1._y,rnode1._z);
//fprintf(E->trace.fpt,"Nodal points: 2: %f %f %f\n",rnode2._x,rnode2._y,rnode2._z);
//fprintf(E->trace.fpt,"Nodal points: 3: %f %f %f\n",rnode3._x,rnode3._y,rnode3._z);
//fprintf(E->trace.fpt,"Nodal points: 4: %f %f %f\n",rnode4._x,rnode4._y,rnode4._z);
fflush(E->trace.fpt);
exit(10);
}
// If any radial component is small, we are on a boundary
if (fabs(rad1)<=tiny||fabs(rad2)<=tiny||fabs(rad3)<=tiny||fabs(rad4)<=tiny)
{
// Convert our test point to spherical coordinates
sc = test_point.toSpherical();
// Ensure it is within bounds
if (sc._theta < 0) sc._theta += 2*M_PI;
sc._rad = 1;
// Nudge the point by eps
if (sc._theta <= M_PI/2.0) sc._theta += eps;
else sc._theta -= eps;
sc._phi += eps;
// Convert the nudged test point back to Cartesian coordinates
test_point = sc.toCartesian();
number_of_tries++;
goto try_again;
}
if (rad1>0.0 && rad2>0.0 && rad3>0.0 && rad4>0.0) return 1;
else return 0;
/*
fprintf(stderr,"%d: icheck: %d\n",E->parallel.me,icheck);
fprintf(stderr,"%d: rads: %f %f %f %f\n",E->parallel.me,rad1,rad2,rad3,rad4);
*/
}
/****************************************************************************/
/* FINDRADIAL */
/* */
/* This function finds the radial component of a Cartesian vector */
static double findradial(CartesianCoord vec,
BoundaryPoint bp)
{
double radialparti,radialpartj,radialpartk;
radialparti=vec._x*bp.sin_theta*bp.cos_phi;
radialpartj=vec._y*bp.sin_theta*bp.sin_phi;
radialpartk=vec._z*bp.cos_theta;
return radialparti+radialpartj+radialpartk;
}
/******************************************************************/
/* FIX ANGLE */
/* */
/* This function constrains the value of angle to be */
/* between 0 and 2 PI */
/* */
static void fix_angle(double *angle)
{
const double two_pi = 2.0*M_PI;
double d2 = floor(*angle / two_pi);
*angle -= two_pi * d2;
}
/********** IGET ELEMENT *****************************************/
/* */
/* This function returns the the real element for a given point. */
/* Returns UNDEFINED_ELEMENT if not in this cap. */
/* Returns -1 if in this cap but cannot find the element. */
/* iprevious_element, if known, is the last known element. If */
/* it is not known, input a negative number. */
int full_iget_element(struct All_variables *E,
int j, int iprevious_element,
CartesianCoord cc,
SphericalCoord sc)
{
int iregel;
int iel;
int ntheta,nphi;
int ival;
int ichoice;
int kk;
int ineighbor;
int icorner[5];
int elx,ely,elz,elxz;
int ifinal_iel;
int nelem;
elx=E->lmesh.elx;
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,sc._rad);
if (ival!=1) return UNDEFINED_ELEMENT;
}
/* do quick search to see if element can be easily found. */
/* note that element may still be out of this cap, but */
/* it is probably fast to do a quick search before */
/* checking cap */
/* get regular element number */
iregel=iget_regel(E,j,sc._theta,sc._phi,&ntheta,&nphi);
if (iregel<=0)
{
return UNDEFINED_ELEMENT;
}
/* AKMA put safety here or in make grid */
if (E->trace.regtoel[j][0][iregel]==0)
{
iel=E->trace.regtoel[j][1][iregel];
goto foundit;
}
/* first check previous element */
if (iprevious_element>0)
{
ival=icheck_element_column(E,j,iprevious_element,cc,sc._rad);
if (ival==1)
{
iel=iprevious_element;
goto foundit;
}
}
/* Check all regular mapping choices */
ichoice=0;
if (E->trace.regtoel[j][0][iregel]>0)
{
ichoice=E->trace.regtoel[j][0][iregel];
for (kk=1;kk<=ichoice;kk++)
{
nelem=E->trace.regtoel[j][kk][iregel];
if (nelem!=iprevious_element)
{
ival=icheck_element_column(E,j,nelem,cc,sc._rad);
if (ival==1)
{
iel=nelem;
goto foundit;
}
}
}
}
/* If here, it means that tracer could not be found quickly with regular element map */
/* First check previous element neighbors */
if (iprevious_element>0)
{
iel=icheck_column_neighbors(E,j,iprevious_element,cc,sc._rad);
if (iel>0)
{
goto foundit;
}
}
/* check if still in cap */
ival=full_icheck_cap(E,0,cc,sc._rad);
if (ival==0) return UNDEFINED_ELEMENT;
/* if here, still in cap (hopefully, without a doubt) */
/* check cap corners (they are sometimes tricky) */
elxz=elx*elz;
icorner[1]=elz;
icorner[2]=elxz;
icorner[3]=elxz*(ely-1)+elz;
icorner[4]=elxz*ely;
for (kk=1;kk<=4;kk++)
{
ival=icheck_element_column(E,j,icorner[kk],cc,sc._rad);
if (ival>0)
{
iel=icorner[kk];
goto foundit;
}
}
/* if previous element is not known, check neighbors of those tried in iquick... */
if (iprevious_element<0)
{
if (ichoice>0)
{
for (kk=1;kk<=ichoice;kk++)
{
ineighbor=E->trace.regtoel[j][kk][iregel];
iel=icheck_column_neighbors(E,j,ineighbor,cc,sc._rad);
if (iel>0)
{
goto foundit;
}
}
}
}
/* As a last resort, check all element columns */
E->trace.istat1++;
iel=icheck_all_columns(E,j,cc,sc._rad);
/*
fprintf(E->trace.fpt,"WARNING(full_iget_element)-doing a full search!\n");
fprintf(E->trace.fpt," Most often means tracers have moved more than 1 element away\n");
fprintf(E->trace.fpt," or regular element resolution is way too low.\n");
fprintf(E->trace.fpt," COLUMN: %d \n",iel);
fprintf(E->trace.fpt," PREVIOUS ELEMENT: %d \n",iprevious_element);
fprintf(E->trace.fpt," x,y,z,theta,phi,rad: %f %f %f %f %f %f\n",x,y,z,theta,phi,rad);
fflush(E->trace.fpt);
if (E->trace.itracer_warnings) exit(10);
*/
if (E->trace.istat1%100==0)
{
fprintf(E->trace.fpt,"Checked all elements %d times already this turn\n",E->trace.istat1);
fflush(E->trace.fpt);
}
if (iel>0)
{
goto foundit;
}
/* if still here, there is a problem */
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",
cc._x,cc._y,cc._z,sc._theta,sc._phi,iregel);
fflush(E->trace.fpt);
return -1;
foundit:
/* find radial element */
ifinal_iel=iget_radial_element(E,j,iel,sc._rad);
return ifinal_iel;
}
/***** IGET RADIAL ELEMENT ***********************************/
/* */
/* This function returns the proper radial element, given */
/* an element (iel) from the column. */
static int iget_radial_element(struct All_variables *E,
int j, int iel,
double rad)
{
int elz=E->lmesh.elz;
int ibottom_element;
int iradial_element;
int node;
int kk;
int idum;
double top_rad;
/* first project to the lowest element in the column */
idum=(iel-1)/elz;
ibottom_element=idum*elz+1;
iradial_element=ibottom_element;
for (kk=1;kk<=elz;kk++)
{
node=E->ien[j][iradial_element].node[8];
top_rad=E->sx[j][3][node];
if (rad<top_rad) goto found_it;
iradial_element++;
} /* end kk */
/* should not be here */
fprintf(E->trace.fpt,"Error(iget_radial_element)-out of range %f %d %d %d\n",rad,j,iel,ibottom_element);
fflush(E->trace.fpt);
exit(10);
found_it:
return iradial_element;
}
/*********** IGET REGEL ******************************************/
/* */
/* This function returns the regular element in which a point */
/* exists. If not found, returns UNDEFINED_ELEMENT. */
/* npi and ntheta are modified for later use */
static int iget_regel(struct All_variables *E, int j,
double theta, double phi,
int *ntheta, int *nphi)
{
int iregel;
int idum;
double rdum;
/* first check whether theta is in range */
if (theta<E->trace.thetamin[j]) return UNDEFINED_ELEMENT;
if (theta>E->trace.thetamax[j]) return UNDEFINED_ELEMENT;
/* get ntheta, nphi on regular mesh */
rdum=theta-E->trace.thetamin[j];
idum=rdum/E->trace.deltheta[j];
*ntheta=idum+1;
rdum=phi-E->trace.phimin[j];
idum=rdum/E->trace.delphi[j];
*nphi=idum+1;
iregel=*ntheta+(*nphi-1)*E->trace.numtheta[j];
/* check range to be sure */
if (iregel>E->trace.numregel[j]) return UNDEFINED_ELEMENT;
if (iregel<1) return UNDEFINED_ELEMENT;
return iregel;
}
/****************************************************************/
/* DEFINE UV SPACE */
/* */
/* This function defines nodal points as orthodrome coordinates */
/* u and v. In uv space, great circles form straight lines. */
/* This is used for interpolation method 1 */
/* E->gnomonic[node].u = u */
/* E->gnomonic[node].v = v */
void define_uv_space(struct All_variables *E)
{
const int j = 1;
const int lev = E->mesh.levmax;
int refnode;
int i, n;
double u, v, cosc, theta_f, phi_f, dphi, cosd;
double *cost, *sint, *cosf, *sinf;
sint = E->SinCos[lev][j][0];
sinf = E->SinCos[lev][j][1];
cost = E->SinCos[lev][j][2];
cosf = E->SinCos[lev][j][3];
/* uv space requires a reference point */
/* use the point at middle of the cap */
refnode = 1 + E->lmesh.noz * ((E->lmesh.noy / 2) * E->lmesh.nox
+ E->lmesh.nox / 2);
phi_f = E->trace.gnomonic_reference_phi = E->sx[j][2][refnode];
/** debug **
theta_f = E->sx[j][1][refnode];
for (i=1; i<=E->lmesh.nsf; i++) {
fprintf(E->trace.fpt, "i=%d (%e %e %e %e)\n",
i, sint[i], sinf[i], cost[i], cosf[i]);
}
fprintf(E->trace.fpt, "%d %d %d ref=(%e %e)\n",
E->lmesh.noz, E->lmesh.nsf, refnode, theta_f, phi_f);
/**/
/* store cos(theta_f) and sin(theta_f) */
E->trace.gnomonic->insert(std::make_pair(0, CoordUV(cost[refnode], sint[refnode])));
/* convert each nodal point to u and v */
for (i=1, n=1; i<=E->lmesh.nsf; i++, n+=E->lmesh.noz) {
dphi = E->sx[j][2][n] - phi_f;
cosd = cos(dphi);
cosc = cost[refnode] * cost[n] + sint[refnode] * sint[n] * cosd;
u = sint[n] * sin(dphi) / cosc;
v = (sint[refnode] * cost[n] - cost[refnode] * sint[n] * cosd)
/ cosc;
E->trace.gnomonic->insert(std::make_pair(i, CoordUV(u, v)));
/** debug **
fprintf(E->trace.fpt, "n=%d ns=%d cosc=%e (%e %e) -> (%e %e)\n",
n, i, cosc, E->sx[j][1][n], E->sx[j][2][n], u, v);
/**/
}
}
/***************************************************************/
/* DETERMINE SHAPE COEFFICIENTS */
/* */
/* An initialization function that determines the coeffiecients*/
/* to all element shape functions. */
/* This method uses standard linear shape functions of */
/* triangular elements. (See Cuvelier, Segal, and */
/* van Steenhoven, 1986). This is all in UV space. */
/* */
/* shape_coefs[cap][wedge][3 shape functions*3 coefs][nelems] */
void determine_shape_coefficients(struct All_variables *E)
{
const int j = 1;
int nelem, iwedge, kk, i;
int snode;
CoordUV uv[5], xy1, xy2, xy3;
/* first, allocate memory */
for(iwedge=0; iwedge<2; iwedge++) {
if ((E->trace.shape_coefs[j][iwedge] =
(TriElemLinearShapeFunc **)malloc((E->lmesh.snel+1)*sizeof(TriElemLinearShapeFunc*))) == NULL) {
fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
fflush(E->trace.fpt);
exit(10);
}
}
for (i=1, nelem=1; i<=E->lmesh.snel; i++, nelem+=E->lmesh.elz) {
/* find u,v of local nodes at one radius */
for(kk=1; kk<=4; kk++) {
snode = (E->ien[j][nelem].node[kk]-1) / E->lmesh.noz + 1;
uv[kk] = (*E->trace.gnomonic)[snode];
}
for(iwedge=0; iwedge<2; iwedge++) {
if (iwedge == 0) {
xy1 = uv[1];
xy2 = uv[2];
xy3 = uv[3];
} else if (iwedge == 1) {
xy1 = uv[1];
xy2 = uv[3];
xy3 = uv[4];
}
E->trace.shape_coefs[j][iwedge][i] = new TriElemLinearShapeFunc(xy1, xy2, xy3);
/** debug **
fprintf(E->trace.fpt, "el=%d els=%d iwedge=%d shape=(%e %e %e, %e %e %e, %e %e %e)\n",
nelem, i, iwedge,
E->trace.shape_coefs[j][iwedge][0][i],
E->trace.shape_coefs[j][iwedge][1][i],
E->trace.shape_coefs[j][iwedge][2][i],
E->trace.shape_coefs[j][iwedge][3][i],
E->trace.shape_coefs[j][iwedge][4][i],
E->trace.shape_coefs[j][iwedge][5][i],
E->trace.shape_coefs[j][iwedge][6][i],
E->trace.shape_coefs[j][iwedge][7][i],
E->trace.shape_coefs[j][iwedge][8][i]);
/**/
} /* end wedge */
} /* end elem */
}
/* &&&&&&&&&&&&&&&&&&&& ANALYTICAL TESTS &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&**************/
/**************** ANALYTICAL TEST *********************************************************/
/* */
/* This function (and the 2 following) are used to test advection of tracers by assigning */
/* a test function (in "analytical_test_function"). */
void analytical_test(struct All_variables *E)
{
#if 0
int kk;
int nsteps;
int j;
int my_number,number;
int nrunge_steps;
int nrunge_refinement;
double dt;
double runge_dt;
double theta,phi,rad;
double time;
double my_theta0,my_phi0,my_rad0;
double my_thetaf,my_phif,my_radf;
double theta0,phi0,rad0;
double thetaf,phif,radf;
SphericalCoord vel_s, x0_s, xf_s;
CartesianCoord vel_c, x0_c, xf_c;
double vec[4];
double runge_path_length,runge_time;
double x0,y0,z0;
double difference;
double difperpath;
TracerList::iterator tr;
fprintf(E->trace.fpt,"Starting Analytical Test\n");
if (E->parallel.me==0) fprintf(stderr,"Starting Analytical Test\n");
fflush(E->trace.fpt);
/* Reset Box cushion to 0 */
E->trace.box_cushion=0.0000;
/* test paramters */
nsteps=200;
dt=0.0001;
E->advection.timestep=dt;
fprintf(E->trace.fpt,"steps: %d dt: %f\n",nsteps,dt);
/* Assign test velocity to Citcom nodes */
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
for (kk=1;kk<=E->lmesh.nno;kk++)
{
theta=E->sx[j][1][kk];
phi=E->sx[j][2][kk];
rad=E->sx[j][3][kk];
analytical_test_function(E,SphericalCoord(theta,phi,rad),vel_s,vel_c);
E->sphere.cap[j].V[1][kk]=vel_s._theta;
E->sphere.cap[j].V[2][kk]=vel_s._phi;
E->sphere.cap[j].V[3][kk]=vel_s._rad;
}
}
time=0.0;
my_theta0=0.0;
my_phi0=0.0;
my_rad0=0.0;
my_thetaf=0.0;
my_phif=0.0;
my_radf=0.0;
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
if (E->trace.tracers[j].size()>10)
{
fprintf(E->trace.fpt,"Warning(analytical)-too many tracers to print!\n");
fflush(E->trace.fpt);
if (E->trace.itracer_warnings) exit(10);
}
}
/* print initial positions */
E->monitor.solution_cycles=0;
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
bool first = true;
for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr)
{
theta=tr->theta();
phi=tr->phi();
rad=tr->rad();
fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
if (first) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
if (first)
{
my_theta0=theta;
my_phi0=phi;
my_rad0=rad;
}
first = false;
}
}
/* advect tracers */
for (kk=1;kk<=nsteps;kk++)
{
E->monitor.solution_cycles=kk;
time += dt;
predict_tracers(E);
correct_tracers(E);
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
bool first = true;
for (tr=E->trace.tracers[j].begin();tr!=E->trace.tracers[j].end();++tr)
{
theta=tr->theta();
phi=tr->phi();
rad=tr->rad();
fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
if (first) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
if ((kk==nsteps) && (first))
{
my_thetaf=theta;
my_phif=phi;
my_radf=rad;
}
}
}
}
/* Get ready for comparison to Runge-Kutte (only works for one tracer) */
fflush(E->trace.fpt);
parallel_process_sync(E);
fprintf(E->trace.fpt,"\n\nComparison to Runge-Kutte\n");
if (E->parallel.me==0) fprintf(stderr,"Comparison to Runge-Kutte\n");
for (j=1;j<=E->sphere.caps_per_proc;j++)
{
my_number=E->trace.tracers[j].size();
}
MPI_Allreduce(&my_number,&number,1,MPI_INT,MPI_SUM,E->parallel.world);
fprintf(E->trace.fpt,"Number of tracers: %d\n", number);
if (E->parallel.me==0) fprintf(stderr,"Number of tracers: %d\n", number);
/* if more than 1 tracer, exit */
if (number!=1)
{
fprintf(E->trace.fpt,"(Note: RK comparison only appropriate for one tracing particle (%d here) \n",number);
if (E->parallel.me==0) fprintf(stderr,"(Note: RK comparison only appropriate for one tracing particle (%d here) \n",number);
fflush(E->trace.fpt);
parallel_process_termination();
}
/* communicate starting and final positions */
MPI_Allreduce(&my_theta0,&theta0,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
MPI_Allreduce(&my_phi0,&phi0,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
MPI_Allreduce(&my_rad0,&rad0,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
MPI_Allreduce(&my_thetaf,&thetaf,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
MPI_Allreduce(&my_phif,&phif,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
MPI_Allreduce(&my_radf,&radf,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
x0_s = SphericalCoord(theta0, phi0, rad0);
nrunge_refinement=1000;
nrunge_steps=nsteps*nrunge_refinement;
runge_dt=dt/(1.0*nrunge_refinement);
analytical_runge_kutte(E,nrunge_steps,runge_dt,x0_s,x0_c,xf_s,xf_c,vec);
runge_time=vec[1];
runge_path_length=vec[2];
/* initial coordinates - both citcom and RK */
CartesianCoord ccf;
SphericalCoord scf(thetaf, phif, radf);
x0=x0_c._x;
y0=x0_c._y;
z0=x0_c._z;
/* convert final citcom coords into cartesian */
ccf = scf.toCartesian();
difference=ccf.dist(xf_c);
difperpath=difference/runge_path_length;
/* Print out results */
fprintf(E->trace.fpt,"Citcom calculation: steps: %d dt: %f\n",nsteps,dt);
fprintf(E->trace.fpt," (nodes per cap: %d x %d x %d)\n",E->lmesh.nox,E->lmesh.noy,(E->lmesh.noz-1)*E->parallel.nprocz+1);
fprintf(E->trace.fpt," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
fprintf(E->trace.fpt," final position: theta: %f phi: %f rad: %f\n", thetaf,phif,radf);
fprintf(E->trace.fpt," (final time: %f) \n",time );
fprintf(E->trace.fpt,"\n\nRunge-Kutte calculation: steps: %d dt: %g\n",nrunge_steps,runge_dt);
fprintf(E->trace.fpt," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
fprintf(E->trace.fpt," final position: theta: %f phi: %f rad: %f\n",xf_s._theta,xf_s._phi,xf_s._rad);
fprintf(E->trace.fpt," path length: %f \n",runge_path_length );
fprintf(E->trace.fpt," (final time: %f) \n",runge_time );
fprintf(E->trace.fpt,"\n\n Difference between Citcom and RK: %e (diff per path length: %e)\n\n",difference,difperpath);
if (E->parallel.me==0)
{
fprintf(stderr,"Citcom calculation: steps: %d dt: %f\n",nsteps,dt);
fprintf(stderr," (nodes per cap: %d x %d x %d)\n",E->lmesh.nox,E->lmesh.noy,(E->lmesh.noz-1)*E->parallel.nprocz+1);
fprintf(stderr," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
fprintf(stderr," final position: theta: %f phi: %f rad: %f\n", thetaf,phif,radf);
fprintf(stderr," (final time: %f) \n",time );
fprintf(stderr,"\n\nRunge-Kutte calculation: steps: %d dt: %f\n",nrunge_steps,runge_dt);
fprintf(stderr," starting position: theta: %f phi: %f rad: %f\n", theta0,phi0,rad0);
fprintf(stderr," final position: theta: %f phi: %f rad: %f\n",xf_s._theta,xf_s._phi,xf_s._rad);
fprintf(stderr," path length: %f \n",runge_path_length );
fprintf(stderr," (final time: %f) \n",runge_time );
fprintf(stderr,"\n\n Difference between Citcom and RK: %e (diff per path length: %e)\n\n",difference,difperpath);
}
fflush(E->trace.fpt);
#endif
}
/*************** ANALYTICAL RUNGE KUTTE ******************/
/* */
void analytical_runge_kutte(struct All_variables *E, int nsteps, double dt,
SphericalCoord &x0_s, CartesianCoord &x0_c,
SphericalCoord &xf_s, CartesianCoord &xf_c, double *vec)
{
int kk;
CartesianCoord cc0, ccp, ccc, vel0_c, velp_c;
SphericalCoord sc0, scp, scc, vel0_s, velp_s;
double time, path;
sc0 = x0_s;
cc0 = sc0.toCartesian();
/* fill initial cartesian vector to send back */
x0_c = cc0;
time=0.0;
path=0.0;
for (kk=0;kk<nsteps;kk++)
{
// Get velocity at initial position
analytical_test_function(E,sc0,vel0_s,vel0_c);
// Find predicted midpoint position
ccp = cc0 + vel0_c * (dt*0.5);
// Convert to spherical
scp = ccp.toSpherical();
// Get velocity at predicted midpoint position
analytical_test_function(E,scp,velp_s,velp_c);
// Find corrected position using midpoint velocity
ccc = cc0 + velp_c * dt;
// Convert to spherical
scc = ccc.toSpherical();
// Compute path length
path += ccc.dist(cc0);
time += dt;
cc0 = ccc;
/* next time step */
}
/* fill final spherical and cartesian vectors to send back */
xf_s = scc;
xf_c = ccc;
vec[1]=time;
vec[2]=path;
}
/**************** ANALYTICAL TEST FUNCTION ******************/
/* */
/* vel_s => velocity in spherical directions */
/* */
/* vel_c => velocity in cartesian directions */
void analytical_test_function(struct All_variables *E, SphericalCoord sc, SphericalCoord &vel_s, CartesianCoord &vel_c)
{
/* This is where the function is given in spherical */
vel_s._theta=50.0*sc._rad*cos(sc._phi);
vel_s._phi=100.0*sc._rad*sin(sc._theta);
vel_s._rad=25.0*sc._rad;
/* Convert the function into cartesian */
vel_c = vel_s.toCartesian();
}
/**** PDEBUG ***********************************************************/
void pdebug(struct All_variables *E, int i)
{
fprintf(E->trace.fpt,"HERE (Before Sync): %d\n",i);
fflush(E->trace.fpt);
parallel_process_sync(E);
fprintf(E->trace.fpt,"HERE (After Sync): %d\n",i);
fflush(E->trace.fpt);
return;
}
