[cig-commits] r6145 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Mar 1 13:02:01 PST 2007
Author: tan2
Date: 2007-03-01 13:02:00 -0800 (Thu, 01 Mar 2007)
New Revision: 6145
Added:
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
Modified:
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Makefile.am
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/tracer_defs.h
Log:
Copied from Allen McNamara's tracer code. The code will fail to compile
Added: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-01 19:36:33 UTC (rev 6144)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-03-01 21:02:00 UTC (rev 6145)
@@ -0,0 +1,7018 @@
+#include<math.h>
+#include<sys/types.h>
+#include "element_definitions.h"
+#include "global_defs.h"
+#include "mpi.h"
+#include<stdlib.h>
+
+#if(0)
+#endif
+
+void pdebug();
+void debug();
+void parallel_process_sync();
+void parallel_process_termination();
+double CPU_time0();
+
+/******* TRACING *************************************************************/
+/* */
+/* This function is the primary tracing routine called from Citcom.c */
+/* In this code, unlike the original 3D cartesian code, force is filled */
+/* during Stokes solution. No need to call thermal_buoyancy() after tracing. */
+
+void tracing(E,itimes_here_this_step)
+ struct All_variables *E;
+ int itimes_here_this_step;
+{
+
+
+static double start_time;
+static double end_time;
+
+void fill_composition();
+void thermal_buoyancy();
+void check_sum();
+void tracer_post_processing();
+void write_tracers();
+void predict_tracers();
+void correct_tracers();
+double CPU_time0();
+void trace_filter();
+
+
+ parallel_process_sync();
+ start_time=CPU_time0();
+
+ if (itimes_here_this_step==1)
+ {
+
+ E->trace.stat_trace_time=0.0;
+
+ fprintf(E->trace.fpt,"STEP %d\n",E->monitor.solution_cycles);
+ fflush(E->trace.fpt);
+
+/* presave last timesteps tracers as preback */
+
+ if ( (E->monitor.solution_cycles % E->control.record_every)==0)
+ {
+ write_tracers(E,3);
+ }
+ }
+ else if (itimes_here_this_step==2)
+ {
+ E->trace.stat_stokes_time=start_time-end_time;
+ }
+
+/* advect tracers */
+
+ if (itimes_here_this_step==1)
+ {
+ if (E->trace.itracer_advection_scheme==1)
+ {
+ predict_tracers(E);
+ correct_tracers(E);
+ }
+ else if (E->trace.itracer_advection_scheme==2)
+ {
+ predict_tracers(E);
+ }
+ }
+ else if (itimes_here_this_step==2)
+ {
+ if (E->trace.itracer_advection_scheme==1)
+ {
+ }
+ else if (E->trace.itracer_advection_scheme==2)
+ {
+ correct_tracers(E);
+ }
+ }
+
+
+
+ if (E->trace.itracer_type==1)
+ {
+ fill_composition(E);
+ }
+
+ end_time=CPU_time0();
+ E->trace.stat_trace_time=E->trace.stat_trace_time+(end_time-start_time);
+
+
+ if (itimes_here_this_step==2)
+ {
+ check_sum(E);
+ if (E->trace.itracer_type==1)
+ {
+ trace_filter(E);
+ }
+ tracer_post_processing(E);
+ }
+
+return;
+}
+
+/********* TRACER POST PROCESSING ****************************************/
+
+void tracer_post_processing(E)
+ struct All_variables *E;
+
+{
+
+char output_file[200];
+
+double convection_time,tracer_time;
+double trace_fraction,total_time;
+
+void get_bulk_composition();
+void write_compositional_field();
+void write_tracers();
+void write_radial_horizontal_averages();
+
+static int been_here=0;
+
+ if (E->trace.itracer_type==1) get_bulk_composition(E);
+
+ if ( ((E->monitor.solution_cycles % E->control.record_every)==0) || (been_here==0))
+ {
+ write_radial_horizontal_averages(E);
+ write_tracers(E,1);
+ if (E->trace.itracer_type==1) write_compositional_field(E);
+ }
+ if ( ((E->monitor.solution_cycles % E->trace.iwrite_tracers_every)==0) &&
+ ((E->monitor.solution_cycles!=E->control.restart)&&(E->monitor.solution_cycles!=0) ))
+ {
+ write_tracers(E,2);
+ }
+
+
+ fprintf(E->trace.fpt,"Number of times for all element search %d\n",E->trace.istat1);
+ if (been_here!=0)
+ {
+ fprintf(E->trace.fpt,"Number of tracers sent to other processors: %d\n",E->trace.istat_isend);
+ fprintf(E->trace.fpt,"Number of times element columns are checked: %d \n",E->trace.istat_elements_checked);
+ if (E->trace.itracer_type==1)
+ {
+ fprintf(E->trace.fpt,"Empty elements filled with old compositional values: %d (%f percent)\n",
+ E->trace.istat_iempty,(100.0*E->trace.istat_iempty/E->lmesh.nel));
+ }
+ }
+
+ fprintf(E->trace.fpt,"Error fraction: %f (composition: %f)\n",E->trace.error_fraction,E->trace.bulk_composition);
+
+
+
+/* reset statistical counters */
+
+ E->trace.istat_isend=0;
+ E->trace.istat_iempty=0;
+ E->trace.istat_elements_checked=0;
+ E->trace.istat1=0;
+
+/* compositional and error fraction data files */
+
+ if (E->parallel.me==0)
+ {
+ if (been_here==0)
+ {
+ if (E->trace.itracer_type==1)
+ {
+ sprintf(output_file,"%s.error_fraction.data",E->control.data_file2);
+ E->trace.fp_error_fraction=fopen(output_file,"w");
+
+ sprintf(output_file,"%s.composition.data",E->control.data_file2);
+ E->trace.fp_composition=fopen(output_file,"w");
+ }
+ }
+
+ if (E->trace.itracer_type==1)
+ {
+ fprintf(E->trace.fp_error_fraction,"%e %e\n",E->monitor.elapsed_time,E->trace.error_fraction);
+ fprintf(E->trace.fp_composition,"%e %e\n",E->monitor.elapsed_time,E->trace.bulk_composition);
+
+ fflush(E->trace.fp_error_fraction);
+ fflush(E->trace.fp_composition);
+ }
+
+ }
+
+
+ convection_time=E->trace.stat_stokes_time;
+ tracer_time=E->trace.stat_trace_time;
+
+ total_time=convection_time+tracer_time;
+ trace_fraction=tracer_time/total_time;
+
+ if (been_here!=0)
+ {
+ fprintf(E->trace.fpt,"CPU time: Convection: %f Tracing: %f (%f percent) Total: %f\n",
+ convection_time,tracer_time,trace_fraction*100.0,total_time);
+
+ }
+
+ fflush(E->trace.fpt);
+
+ if (been_here==0) been_here++;
+
+return;
+}
+
+/*********** WRITE RADIAL HORIZONTAL AVERAGES ***************************/
+/* */
+/* This function writes the radial distribution of tracers */
+/* (horizontally averaged) */
+
+void write_radial_horizontal_averages(E)
+ struct All_variables *E;
+{
+
+char output_file[200];
+
+int j;
+int kk;
+double halfpoint;
+double *reltrac[13];
+
+static int been_here=0;
+
+void return_horiz_ave();
+void return_elementwise_horiz_ave();
+
+FILE *fp2;
+
+ if (been_here==0)
+ {
+ E->trace.Have_C=(double *)malloc((E->lmesh.noz+2)*sizeof(double));
+ E->trace.Havel_tracers=(double *)malloc((E->lmesh.elz+2)*sizeof(double));
+ }
+
+/* Tracers */
+
+/* first, change from int to double */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ reltrac[j]=(double *) malloc((E->lmesh.nel+1)*sizeof(double));
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ reltrac[j][kk]=(1.0*E->trace.ieltrac[j][kk]);
+ }
+ }
+
+ return_elementwise_horiz_ave(E,reltrac,E->trace.Havel_tracers);
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ free(reltrac[j]);
+ }
+
+ if (E->parallel.me<E->parallel.nprocz)
+ {
+ sprintf(output_file,"%s.ave_tracers.%d.%d",E->control.data_file2,E->parallel.me,E->monitor.solution_cycles);
+ fp2=fopen(output_file,"w");
+ for(kk=1;kk<=E->lmesh.elz;kk++)
+ {
+ halfpoint=0.5*(E->sx[1][3][kk+1]+E->sx[1][3][kk]);
+ fprintf(fp2,"%.4e %.4e\n",halfpoint,E->trace.Havel_tracers[kk]);
+ }
+ fclose(fp2);
+ }
+
+/* Composition */
+
+ if (E->trace.itracer_type==1)
+ {
+ return_horiz_ave(E,E->trace.comp_node,E->trace.Have_C);
+
+
+ if (E->parallel.me<E->parallel.nprocz)
+ {
+ sprintf(output_file,"%s.ave_c.%d.%d",E->control.data_file2,E->parallel.me,E->monitor.solution_cycles);
+ fp2=fopen(output_file,"w");
+ for(kk=1;kk<=E->lmesh.noz;kk++)
+ {
+ fprintf(fp2,"%.4e %.4e\n",E->sx[1][3][kk],E->trace.Have_C[kk]);
+ }
+ fclose(fp2);
+
+ }
+ }
+
+been_here++;
+
+return;
+}
+
+/*********** WRITE TRACERS **********************************************/
+/* */
+/* if iflag=1, rewrite over last save array */
+/* if iflag=2, write periodic save file and gzip */
+
+
+
+void write_tracers(E,iflag)
+ struct All_variables *E;
+ int iflag;
+
+{
+
+char output_file[200];
+/* char doit_string[200]; */
+FILE *fpcomp;
+
+int kk,j;
+double comp;
+double theta,phi,rad;
+
+ if (iflag==1)
+ {
+ sprintf(output_file,"%s.tracers.%d",E->control.data_file2,
+ E->parallel.me);
+ }
+ if (iflag==2)
+ {
+ sprintf(output_file,"%s.tracers.%d.%d",E->control.data_file2,
+ E->parallel.me, E->monitor.solution_cycles);
+ }
+ if (iflag==3)
+ {
+ sprintf(output_file,"%s.tracers.%d.preback",E->control.data_file2,
+ E->parallel.me);
+ }
+ fpcomp=fopen(output_file,"w");
+
+/* This may be slightly incompatible with previous version */
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ fprintf(fpcomp,"%6d %6d %.5e %.5e %.5e %d\n",E->trace.itrac[j][0],
+ E->monitor.solution_cycles,E->monitor.elapsed_time,
+ E->trace.bulk_composition,
+ E->trace.initial_bulk_composition,j);
+
+ comp=0.0;
+ for (kk=1;kk<=E->trace.itrac[j][0];kk++)
+ {
+ theta=E->trace.rtrac[j][0][kk];
+ phi=E->trace.rtrac[j][1][kk];
+ rad=E->trace.rtrac[j][2][kk];
+
+ if (E->trace.itracer_type==1)
+ {
+ comp=E->trace.etrac[j][0][kk];
+ fprintf(fpcomp,"%.12e %.12e %.12e %.12e \n",theta,phi,rad,comp);
+ fflush(fpcomp);
+ }
+ else if (E->trace.itracer_type==0)
+ {
+ fprintf(fpcomp,"%.12e %.12e %.12e \n",theta,phi,rad);
+ fflush(fpcomp);
+ }
+ }
+ fclose(fpcomp);
+
+/* maybe some problem with zipping on the fly? */
+/*
+ if (iflag==2)
+ {
+ sprintf(doit_string,"gzip %s",output_file);
+ system(doit_string);
+ }
+*/
+
+ }
+
+return;
+}
+
+
+
+
+
+/********** WRITE COMPOSITIONAL FIELD ***********************************/
+
+void write_compositional_field(E)
+ struct All_variables *E;
+{
+
+char output_file1[1000];
+char output_file2[1000];
+
+int j,i;
+
+FILE *fp1;
+FILE *fp2;
+
+
+/* comp */
+
+ sprintf(output_file1,"%s.comp.%d.%d",E->control.data_file2,E->parallel.me,E->monitor.solution_cycles);
+
+ fp1=fopen(output_file1,"w");
+
+ fprintf(fp1,"%d %d %.5e %.5e %.5e\n",E->monitor.solution_cycles,E->lmesh.nno,E->monitor.elapsed_time,
+ E->trace.initial_bulk_composition,E->trace.bulk_composition);
+
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ fprintf(fp1,"%3d %7d\n",j,E->lmesh.nno);
+
+ for(i=1;i<=E->lmesh.nno;i++)
+ {
+ fprintf(fp1,"%.6e\n",E->trace.comp_node[j][i]);
+ }
+ }
+
+ fflush(fp1);
+
+/* OLDCOMP */
+
+ sprintf(output_file2,"%s.OLDCOMP.%d",E->control.data_file2,E->parallel.me);
+ fp2=fopen(output_file2,"w");
+ fprintf(fp2,"%d %d %.5e \n",E->lmesh.nel,E->monitor.solution_cycles,E->monitor.elapsed_time);
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ fprintf(fp2,"%d %d\n",j,E->lmesh.nel);
+ for(i=1;i<=E->lmesh.nel;i++)
+ {
+ fprintf(fp2,"%.6e\n",E->trace.oldel[j][i]);
+ }
+ }
+ fflush(fp2);
+
+ if ( ((E->monitor.solution_cycles % E->trace.iwrite_tracers_every)==0) &&
+ ((E->monitor.solution_cycles!=E->control.restart)&&(E->monitor.solution_cycles!=0) ))
+ {
+ sprintf(output_file2,"%s.OLDCOMP.%d.%d",E->control.data_file2,E->parallel.me,E->monitor.solution_cycles);
+ fp2=fopen(output_file2,"w");
+ fprintf(fp2,"%d %d %.5e \n",E->lmesh.nel,E->monitor.solution_cycles,E->monitor.elapsed_time);
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ fprintf(fp2,"%d %d\n",j,E->lmesh.nel);
+ for(i=1;i<=E->lmesh.nel;i++)
+ {
+ fprintf(fp2,"%.6e\n",E->trace.oldel[j][i]);
+ }
+ }
+
+ }
+ fflush(fp2);
+
+
+
+return;
+}
+
+/*********** GET BULK COMPOSITION *******************************/
+
+void get_bulk_composition(E)
+ struct All_variables *E;
+
+{
+
+char output_file[200];
+char input_s[1000];
+
+double return_bulk_value_d();
+double volume;
+double rdum1;
+double rdum2;
+double rdum3;
+
+int ival=0;
+int idum1;
+int istep;
+
+
+FILE *fp;
+
+static int been_here=0;
+
+
+/* ival=0 returns integral not average */
+
+ volume=return_bulk_value_d(E,E->trace.comp_node,ival);
+
+ E->trace.bulk_composition=volume;
+
+/* Here we assume if restart = -1 or 0 tracers are reset */
+/* if restart = 1 tracers may or may not be reset */
+/* (read initial composition from file) */
+
+
+ if (been_here==0)
+ {
+ if (E->trace.ireset_initial_composition==1)
+ {
+ E->trace.initial_bulk_composition=volume;
+ }
+ else
+ {
+
+ if (E->trace.itracer_restart!=1)
+ {
+ fprintf(E->trace.fpt,"ERROR(bulk composition)-wrong reset,restart combo\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,
+ E->parallel.me,E->monitor.solution_cycles);
+
+ fp=fopen(output_file,"r");
+ fgets(input_s,200,fp);
+ sscanf(input_s,"%d %d %lf %lf %lf",
+ &istep,&idum1,&rdum1,&rdum2,&rdum3);
+
+ E->trace.initial_bulk_composition=rdum2;
+ fclose(fp);
+
+ if (istep!=E->monitor.solution_cycles)
+ {
+ fprintf(E->trace.fpt,"ERROR(get_bulk_composition) %d %d\n",
+ istep,E->monitor.solution_cycles);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+ E->trace.error_fraction=((volume-E->trace.initial_bulk_composition)/
+ E->trace.initial_bulk_composition);
+
+ parallel_process_sync();
+
+been_here++;
+return;
+}
+
+
+/*********** PREDICT TRACERS **********************************************/
+/* */
+/* This function predicts tracers performing an euler step */
+/* */
+/* */
+/* Note positions used in tracer array */
+/* [positions 0-5 are always fixed with current coordinates */
+/* regardless of advection scheme]. */
+/* Positions 6-8 contain original Cartesian coordinates. */
+/* Positions 9-11 contain original Cartesian velocities. */
+/* */
+
+
+void predict_tracers(E)
+ struct All_variables *E;
+{
+
+int numtracers;
+int j;
+int kk;
+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];
+
+void get_velocity();
+void get_cartesian_velocity_field();
+void cart_to_sphere();
+void keep_in_sphere();
+void find_tracers();
+
+static int been_here=0;
+
+ dt=E->advection.timestep;
+
+/* if advection scheme is 2, don't have to calculate cartesian velocity again */
+/* (already did after last stokes calculation, unless this is first step) */
+
+ if ((been_here==0) && (E->trace.itracer_advection_scheme==2))
+ {
+ get_cartesian_velocity_field(E);
+ been_here++;
+ }
+
+ if (E->trace.itracer_advection_scheme==1) get_cartesian_velocity_field(E);
+
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+ numtracers=E->trace.itrac[j][0];
+
+ for (kk=1;kk<=numtracers;kk++)
+ {
+
+ theta0=E->trace.rtrac[j][0][kk];
+ phi0=E->trace.rtrac[j][1][kk];
+ rad0=E->trace.rtrac[j][2][kk];
+ x0=E->trace.rtrac[j][3][kk];
+ y0=E->trace.rtrac[j][4][kk];
+ z0=E->trace.rtrac[j][5][kk];
+
+ nelem=E->trace.itrac[j][kk];
+ 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);
+ keep_in_sphere(E,&x_pred,&y_pred,&z_pred,&theta_pred,&phi_pred,&rad_pred);
+
+/* Current Coordinates are always kept in positions 0-5. */
+
+ E->trace.rtrac[j][0][kk]=theta_pred;
+ E->trace.rtrac[j][1][kk]=phi_pred;
+ E->trace.rtrac[j][2][kk]=rad_pred;
+ E->trace.rtrac[j][3][kk]=x_pred;
+ E->trace.rtrac[j][4][kk]=y_pred;
+ E->trace.rtrac[j][5][kk]=z_pred;
+
+/* Fill in original coords (positions 6-8) */
+
+ E->trace.rtrac[j][6][kk]=x0;
+ E->trace.rtrac[j][7][kk]=y0;
+ E->trace.rtrac[j][8][kk]=z0;
+
+/* Fill in original velocities (positions 9-11) */
+
+ E->trace.rtrac[j][9][kk]=velocity_vector[1]; /* Vx */
+ E->trace.rtrac[j][10][kk]=velocity_vector[2]; /* Vy */
+ E->trace.rtrac[j][11][kk]=velocity_vector[3]; /* Vz */
+
+
+ } /* end kk, 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 */
+/* regardless of advection scheme]. */
+/* Positions 6-8 contain original Cartesian coordinates. */
+/* Positions 9-11 contain original Cartesian velocities. */
+/* */
+
+
+void correct_tracers(E)
+ struct All_variables *E;
+{
+
+int j;
+int kk;
+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;
+
+void get_velocity();
+void get_cartesian_velocity_field();
+void cart_to_sphere();
+void keep_in_sphere();
+void find_tracers();
+
+
+ dt=E->advection.timestep;
+
+ if ((E->parallel.me==0) && (E->trace.itracer_advection_scheme==2) )
+ {
+ fprintf(stderr,"AA:Correcting Tracers\n");
+ fflush(stderr);
+ }
+
+
+/* If scheme==1, use same velocity (t=0) */
+/* Else if scheme==2, use new velocity (t=dt) */
+
+ if (E->trace.itracer_advection_scheme==2)
+ {
+ get_cartesian_velocity_field(E);
+ }
+
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+ for (kk=1;kk<=E->trace.itrac[j][0];kk++)
+ {
+
+ theta_pred=E->trace.rtrac[j][0][kk];
+ phi_pred=E->trace.rtrac[j][1][kk];
+ rad_pred=E->trace.rtrac[j][2][kk];
+ x_pred=E->trace.rtrac[j][3][kk];
+ y_pred=E->trace.rtrac[j][4][kk];
+ z_pred=E->trace.rtrac[j][5][kk];
+
+ x0=E->trace.rtrac[j][6][kk];
+ y0=E->trace.rtrac[j][7][kk];
+ z0=E->trace.rtrac[j][8][kk];
+
+ Vx0=E->trace.rtrac[j][9][kk];
+ Vy0=E->trace.rtrac[j][10][kk];
+ Vz0=E->trace.rtrac[j][11][kk];
+
+ nelem=E->trace.itrac[j][kk];
+
+ 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);
+ keep_in_sphere(E,&x_cor,&y_cor,&z_cor,&theta_cor,&phi_cor,&rad_cor);
+
+/* Fill in Current Positions (other positions are no longer important) */
+
+ E->trace.rtrac[j][0][kk]=theta_cor;
+ E->trace.rtrac[j][1][kk]=phi_cor;
+ E->trace.rtrac[j][2][kk]=rad_cor;
+ E->trace.rtrac[j][3][kk]=x_cor;
+ E->trace.rtrac[j][4][kk]=y_cor;
+ E->trace.rtrac[j][5][kk]=z_cor;
+
+ } /* end kk, correcting tracers */
+} /* end caps */
+/* find new tracer elements and caps */
+
+ find_tracers(E);
+
+return;
+}
+
+/******** GET VELOCITY ***************************************/
+
+void get_velocity(E,j,nelem,theta,phi,rad,velocity_vector)
+ struct All_variables *E;
+ int j,nelem;
+ double theta,phi,rad;
+ double *velocity_vector;
+{
+
+void gnomonic_interpolation();
+
+/* gnomonic projection */
+
+ if (E->trace.itracer_interpolation_scheme==1)
+ {
+ gnomonic_interpolation(E,j,nelem,theta,phi,rad,velocity_vector);
+ }
+ else if (E->trace.itracer_interpolation_scheme==2)
+ {
+ fprintf(E->trace.fpt,"Error(get velocity)-not ready for simple average interpolation scheme\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ else
+ {
+ fprintf(E->trace.fpt,"Error(get velocity)-not ready for other interpolation schemes\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+return;
+}
+
+/********************** GNOMONIC INTERPOLATION *********************************/
+/* */
+/* This function interpolates tracer velocity using gnominic interpolation. */
+/* 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. */
+/* The transformations require a reference point along each cap. Here, the */
+/* midpoint is used. */
+/* 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 gnomonic_interpolation(E,j,nelem,theta,phi,rad,velocity_vector)
+ struct All_variables *E;
+ int j,nelem;
+ double theta,phi,rad;
+ double *velocity_vector;
+{
+
+int iwedge,inum;
+int kk;
+int ival;
+int itry;
+
+double u,v;
+double shape2d[4];
+double shaperad[3];
+double shape[7];
+double vx[7],vy[7],vz[7];
+double x,y,z;
+
+
+int maxlevel=E->mesh.levmax;
+
+static double eps=-1e-4;
+
+void get_radial_shape();
+void sphere_to_cart();
+void spherical_to_uv();
+void get_2dshape();
+int iget_element();
+int icheck_element();
+int icheck_column_neighbors();
+
+
+/* find u and v using spherical coordinates */
+
+ spherical_to_uv(E,j,theta,phi,&u,&v);
+
+ inum=0;
+ itry=1;
+
+try_again:
+
+/* Check first wedge (1 of 2) */
+
+ iwedge=1;
+
+next_wedge:
+
+/* determine shape functions of wedge */
+/* There are 3 shape functions for the triangular wedge */
+
+ get_2dshape(E,j,nelem,u,v,iwedge,shape2d);
+
+/* if any shape functions are negative, goto next wedge */
+
+ if (shape2d[1]<eps||shape2d[2]<eps||shape2d[3]<eps)
+ {
+ inum=inum+1;
+/* 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)
+ {
+ fprintf(E->trace.fpt,"ERROR(gnomonic_interpolation)-inum>1\n");
+ fprintf(E->trace.fpt,"shape %f %f %f\n",shape2d[1],shape2d[2],shape2d[3]);
+ fprintf(E->trace.fpt,"u %f v %f element: %d \n",u,v, nelem);
+ fprintf(E->trace.fpt,"Element uv boundaries: \n");
+ for(kk=1;kk<=4;kk++)
+ fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,E->trace.UV[j][1][E->ien[j][nelem].node[kk]],E->trace.UV[j][2][E->ien[j][nelem].node[kk]]);
+ fprintf(E->trace.fpt,"theta: %f phi: %f rad: %f\n",theta,phi,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]]);
+ 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=iget_element(E,j,-99,x,y,z,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);
+ nelem=ival;
+ ival=icheck_element(E,j,nelem,x,y,z,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=2;
+ goto next_wedge;
+ }
+
+/* Determine radial shape functions */
+/* There are 2 shape functions radially */
+
+ get_radial_shape(E,j,nelem,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[1]=shaperad[1]*shape2d[1];
+ shape[2]=shaperad[1]*shape2d[2];
+ shape[3]=shaperad[1]*shape2d[3];
+ shape[4]=shaperad[2]*shape2d[1];
+ shape[5]=shaperad[2]*shape2d[2];
+ shape[6]=shaperad[2]*shape2d[3];
+
+/* depending on wedge, set up velocity points */
+
+ if (iwedge==1)
+ {
+ vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
+ vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
+ vx[2]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[2]];
+ vx[3]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[3]];
+ vx[4]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[5]];
+ vx[5]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[6]];
+ vx[6]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[7]];
+ vy[1]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[1]];
+ vy[2]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[2]];
+ vy[3]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[3]];
+ vy[4]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[5]];
+ vy[5]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[6]];
+ vy[6]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[7]];
+ vz[1]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[1]];
+ vz[2]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[2]];
+ vz[3]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[3]];
+ vz[4]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[5]];
+ vz[5]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[6]];
+ vz[6]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[7]];
+ }
+ if (iwedge==2)
+ {
+ vx[1]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[1]];
+ vx[2]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[3]];
+ vx[3]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[4]];
+ vx[4]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[5]];
+ vx[5]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[7]];
+ vx[6]=E->trace.V0_cart[j][1][E->IEN[maxlevel][j][nelem].node[8]];
+ vy[1]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[1]];
+ vy[2]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[3]];
+ vy[3]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[4]];
+ vy[4]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[5]];
+ vy[5]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[7]];
+ vy[6]=E->trace.V0_cart[j][2][E->IEN[maxlevel][j][nelem].node[8]];
+ vz[1]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[1]];
+ vz[2]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[3]];
+ vz[3]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[4]];
+ vz[4]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[5]];
+ vz[5]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[7]];
+ vz[6]=E->trace.V0_cart[j][3][E->IEN[maxlevel][j][nelem].node[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;
+}
+
+/***************************************************************/
+/* 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). */
+
+
+void get_2dshape(E,j,nelem,u,v,iwedge,shape2d)
+ struct All_variables *E;
+ int j,nelem,iwedge;
+ double u,v;
+ double * shape2d;
+
+{
+
+ double a0,a1,a2;
+
+/* shape function 1 */
+
+ a0=E->trace.shape_coefs[j][iwedge][1][nelem];
+ a1=E->trace.shape_coefs[j][iwedge][2][nelem];
+ a2=E->trace.shape_coefs[j][iwedge][3][nelem];
+
+ shape2d[1]=a0+a1*u+a2*v;
+
+/* shape function 2 */
+
+ a0=E->trace.shape_coefs[j][iwedge][4][nelem];
+ a1=E->trace.shape_coefs[j][iwedge][5][nelem];
+ a2=E->trace.shape_coefs[j][iwedge][6][nelem];
+
+ shape2d[2]=a0+a1*u+a2*v;
+
+/* shape function 3 */
+
+ a0=E->trace.shape_coefs[j][iwedge][7][nelem];
+ a1=E->trace.shape_coefs[j][iwedge][8][nelem];
+ a2=E->trace.shape_coefs[j][iwedge][9][nelem];
+
+ shape2d[3]=a0+a1*u+a2*v;
+
+return;
+}
+
+/***************************************************************/
+/* GET RADIAL SHAPE */
+/* */
+/* This function determines radial shape functions at rad */
+
+void get_radial_shape(E,j,nelem,rad,shaperad)
+ struct All_variables *E;
+ int j,nelem;
+ double rad;
+ double * shaperad;
+
+{
+
+ int node1,node5;
+ double rad1,rad5,f1,f2,delrad;
+
+ static double eps=1e-6;
+ double top_bound=1.0+eps;
+ double bottom_bound=0.0-eps;
+
+ node1=E->IEN[E->mesh.levmax][j][nelem].node[1];
+ node5=E->IEN[E->mesh.levmax][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[1]=1.0-f1=1.0-(1.0-f2)=f2;
+ shaperad[2]=1.0-f2=1.0-(10-f1)=f1;
+*/
+
+ shaperad[1]=f2;
+ shaperad[2]=f1;
+
+/* Some error control */
+
+ if (shaperad[1]>(top_bound)||shaperad[1]<(bottom_bound)||
+ shaperad[2]>(top_bound)||shaperad[2]<(bottom_bound))
+ {
+ fprintf(E->trace.fpt,"ERROR(get_radial_shape)\n");
+ fprintf(E->trace.fpt,"shaperad[1]: %f \n",shaperad[1]);
+ fprintf(E->trace.fpt,"shaperad[2]: %f \n",shaperad[2]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+return;
+}
+
+
+
+
+
+/**************************************************************/
+/* SPHERICAL TO UV */
+/* */
+/* This function transforms theta and phi to new coords */
+/* u and v using gnomonic projection. */
+
+void spherical_to_uv(E,j,theta,phi,u,v)
+ struct All_variables *E;
+ int j;
+ double theta,phi;
+ double *u;
+ double *v;
+
+{
+
+ double theta_f;
+ 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 at the midpoint of the cap */
+
+ theta_f=E->trace.UV[j][1][0];
+ phi_f=E->trace.UV[j][2][0];
+
+ cos_theta_f=E->trace.cos_theta_f;
+ sin_theta_f=E->trace.sin_theta_f;
+
+ cost=cos(theta);
+/*
+ sint=sin(theta);
+*/
+ sint=sqrt(1.0-cost*cost);
+
+ cosp2=cos(phi-phi_f);
+ sinp2=sin(phi-phi_f);
+
+ cosc=cos_theta_f*cost+sin_theta_f*sint*cosp2;
+ cosc=1.0/cosc;
+
+ *u=sint*sinp2*cosc;
+ *v=(sin_theta_f*cost-cos_theta_f*sint*cosp2)*cosc;
+
+return;
+}
+
+
+
+
+/******* TRACE INSTRUCTIONS *********************/
+
+void trace_instructions(E)
+ struct All_variables *E;
+
+{
+
+
+int input_int();
+int input_string();
+int input_double();
+
+/* itracer_type=0 passive */
+/* itracer_type=1 active */
+
+ E->trace.itracer_type=1;
+ input_int("tracer_type",&(E->trace.itracer_type),"1,0,nomax");
+
+/* Restart options */
+
+/* itracer_restart=-1 (read from file) */
+/* itracer_restart=0 (generate array) */
+/* itracer_restart=1 (read from scratch disks) */
+
+ E->trace.itracer_restart=0;
+ input_int("tracer_restart",&(E->trace.itracer_restart),"0,0,nomax");
+
+/* Active tracer inputs */
+
+ if (E->trace.itracer_type==1)
+ {
+
+ E->trace.buoyancy_ratio=1.0;
+ input_double("buoyancy_ratio",&(E->trace.buoyancy_ratio),"1.0");
+
+ E->trace.ireset_initial_composition=1;
+ if (E->trace.itracer_restart==1)
+ {
+ E->trace.ireset_initial_composition=0;
+ input_int("reset_initial_composition",&(E->trace.ireset_initial_composition),"0");
+ }
+
+/* compositional rheology */
+
+/* ibuoy_type=1 (ratio method) */
+
+ E->trace.ibuoy_type=1;
+ input_int("buoy_type",&(E->trace.ibuoy_type),"1,0,nomax");
+ if (E->trace.ibuoy_type!=1)
+ {
+ fprintf(stderr,"Terror-Sorry, only ratio method allowed now\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
+
+/* filtering parameters */
+
+ E->trace.itrace_filter=0;
+ E->trace.itrace_filter_every=1;
+ E->trace.trace_filter_uplimit=0.9;
+ E->trace.trace_filter_downlimit=0.1;
+
+ input_int("trace_filter",&(E->trace.itrace_filter),"1,0,nomax");
+ input_int("trace_filter_every",&(E->trace.itrace_filter_every),"1,0,nomax");
+ if (E->trace.itrace_filter == 1)
+ {
+ input_double("trace_filter_uplimit",&(E->trace.trace_filter_uplimit),"1.0");
+ input_double("trace_filter_downlimit",&(E->trace.trace_filter_downlimit),"1.0");
+ }
+
+ }
+
+/* icompositional_rheology=0 (off) */
+/* icompositional_rheology=1 (on) */
+
+ E->trace.icompositional_rheology=0;
+ input_int("compositional_rheology",
+ &(E->trace.icompositional_rheology),"1,0,nomax");
+
+ E->trace.compositional_rheology_prefactor=1.0;
+
+ if (E->trace.icompositional_rheology>0)
+ {
+ input_double("compositional_prefactor",
+ &(E->trace.compositional_rheology_prefactor),"1.0");
+ }
+
+/* ibuoy_type=0 (absolute method) */
+
+/* icartesian_or_spherical=0 (cartesian coordinate input) */
+/* icartesian_or_spherical=1 (spherical coordinate input) */
+
+ E->trace.icartesian_or_spherical_input=0;
+ if (E->trace.itracer_restart==-1)
+ {
+ sprintf(E->trace.tracer_file,"tracer.in");
+ input_string("tracer_file",E->trace.tracer_file,"tracer.in");
+
+ input_int("cartesian_or_spherical_input",&(E->trace.icartesian_or_spherical_input),"10,0,nomax");
+
+ }
+
+
+ if (E->trace.itracer_restart==0)
+ {
+ E->trace.itperel=10;
+ input_int("tracers_per_element",&(E->trace.itperel),"10,0,nomax");
+
+ E->trace.z_interface=0.5;
+ input_double("z_interface",&(E->trace.z_interface),"0.5");
+ }
+
+ if (E->trace.itracer_restart==1)
+ {
+ }
+
+/* Advection Scheme */
+
+/* itracer_advection_scheme=1 (simple predictor corrector -uses only V(to)) */
+/* itracer_advection_scheme=2 (predictor-corrector - uses V(to) and V(to+dt)) */
+
+ E->trace.itracer_advection_scheme=2;
+ input_int("tracer_advection_scheme",&(E->trace.itracer_advection_scheme),
+ "100,0,nomax");
+
+ if (E->trace.itracer_advection_scheme==1)
+ {
+ }
+ else if (E->trace.itracer_advection_scheme==2)
+ {
+ }
+ else
+ {
+ fprintf(stderr,"Sorry, only advection scheme 1 and 2 available (%d)\n",E->trace.itracer_advection_scheme);
+ fflush(stderr);
+ parallel_process_termination();
+ }
+
+
+/* Output Options */
+
+ E->trace.iwrite_tracers_every=1000000;
+ input_int("write_tracers_every",&(E->trace.iwrite_tracers_every),"1000000,0,nomax");
+
+/* Interpolation Scheme */
+/* itracer_interpolation_scheme=1 (gnometric projection) */
+/* itracer_interpolation_scheme=2 (simple average) */
+
+ E->trace.itracer_interpolation_scheme=1;
+ input_int("tracer_interpolation_scheme",&(E->trace.itracer_interpolation_scheme),
+ "100,0,nomax");
+ if (E->trace.itracer_interpolation_scheme<1 || E->trace.itracer_interpolation_scheme>2)
+ {
+ fprintf(stderr,"Sorry, only interpolation scheme 1 and 2 available\n");
+ fflush(stderr);
+ parallel_process_termination();
+ }
+
+/* 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");
+ input_double("regular_grid_delphi",&(E->trace.delphi[0]),"1.0");
+
+
+/* Analytical Test Function */
+
+ E->trace.ianalytical_tracer_test=0;
+ input_int("analytical_tracer_test",&(E->trace.ianalytical_tracer_test),
+ "100,0,nomax");
+
+return;
+}
+
+/***** INITIALIZE TRACE ************************/
+
+void initialize_trace(E)
+ struct All_variables *E;
+{
+
+char output_file[200];
+
+void write_trace_instructions();
+void viscosity_checks();
+void make_tracer_array();
+void setup_shared_cap_information();
+void initialize_old_composition();
+void fill_composition();
+void find_tracers();
+void setup_dsincos();
+void make_regular_grid();
+void initialize_tracer_elements();
+void define_uv_space();
+void determine_shape_coefficients();
+void check_sum();
+void read_tracer_file();
+void analytical_test();
+void tracer_post_processing();
+void restart_tracers();
+
+
+/* some obscure initial parameters */
+
+ E->trace.box_cushion=0.00001;
+
+/* AKMA turn this back on after debugging */
+ E->trace.itracer_warnings=1;
+
+/* Determine number of tracer quantities */
+
+/* advection_quantites - those needed for advection */
+
+ if (E->trace.itracer_advection_scheme==1) E->trace.number_of_advection_quantities=12;
+ if (E->trace.itracer_advection_scheme==2) E->trace.number_of_advection_quantities=12;
+
+/* extra_quantities - used for composition, etc. */
+/* (can be increased for additional science i.e. tracing chemistry */
+
+ if (E->trace.itracer_type==1) E->trace.number_of_extra_quantities=1;
+ else E->trace.number_of_extra_quantities=0;
+
+ E->trace.number_of_tracer_quantities=E->trace.number_of_advection_quantities +
+ E->trace.number_of_extra_quantities;
+
+
+
+/* Fixed positions in tracer array */
+/* Comp is always in etrac position 0 */
+/* Current coordinates are always kept in rtrac positions 0-5 */
+/* Other positions may be used depending on advection scheme and/or science being done */
+
+/* open tracing output file */
+
+ sprintf(output_file,"%s.tracer_log.%d",E->control.data_file2,
+ E->parallel.me);
+ E->trace.fpt=fopen(output_file,"w");
+
+/* reset time counters */
+
+ E->trace.stat_trace_time=0.0;
+ E->trace.stat_stokes_time=0.0;
+
+
+/* Some error control regarding size of pointer arrays */
+
+ if (E->trace.number_of_advection_quantities>99)
+ {
+ fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rtrac in global.h\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ if (E->trace.number_of_extra_quantities>99)
+ {
+ fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of etrac in global.h\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ if (E->trace.number_of_tracer_quantities>99)
+ {
+ fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in global.h\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+ write_trace_instructions(E);
+
+/* Some error control */
+
+ if (E->sphere.caps_per_proc>1)
+ {
+ fprintf(E->trace.fpt,"This code has not been tested or tried for multiple caps per processor!\n");
+ fprintf(E->trace.fpt,"I would certainly do some testing first \n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+ if (E->trace.icompositional_rheology>0) viscosity_checks(E);
+
+ setup_shared_cap_information(E);
+ if (E->trace.itracer_restart==0) make_tracer_array(E);
+ else if (E->trace.itracer_restart==-1) read_tracer_file(E);
+ else if (E->trace.itracer_restart==1) restart_tracers(E);
+ else
+ {
+ fprintf(E->trace.fpt,"Not ready for other inputs yet\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+/* flush and wait for not real reason but it can't hurt */
+ fflush(E->trace.fpt);
+ parallel_process_sync();
+
+
+ setup_dsincos(E);
+
+ if (E->trace.itracer_interpolation_scheme==1)
+ {
+ define_uv_space(E);
+ determine_shape_coefficients(E);
+ }
+
+
+/* flush and wait for not real reason but it can't hurt */
+ fflush(E->trace.fpt);
+ parallel_process_sync();
+
+ make_regular_grid(E);
+
+/* flush and wait for not real reason but it can't hurt */
+ fflush(E->trace.fpt);
+ parallel_process_sync();
+
+
+/* find elements */
+
+ find_tracers(E);
+
+ check_sum(E);
+
+ if (E->trace.itracer_type==1 && E->trace.ibuoy_type==1)
+ {
+ initialize_old_composition(E);
+ fill_composition(E);
+ }
+
+ if (E->trace.ianalytical_tracer_test==1)
+ {
+ analytical_test(E);
+ parallel_process_termination();
+ }
+
+ tracer_post_processing(E);
+
+return;
+}
+
+/**************** INITIALIZE TRACER ARRAYS ************************************/
+/* */
+/* This function allocates memories to tracer arrays. */
+
+void initialize_tracer_arrays(E,j,number_of_tracers)
+ struct All_variables *E;
+ int j, number_of_tracers;
+{
+
+int kk;
+
+/* itracsize is physical size of tracer array */
+/* (initially make it 25% larger than required */
+
+ E->trace.itracsize[j]=number_of_tracers+number_of_tracers/4;
+
+/* make tracer arrays */
+
+ if ((E->trace.itrac[j]=(int *) malloc(E->trace.itracsize[j]*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(make tracer array)-no memory 1a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ E->trace.itrac[j][0]=0;
+
+ for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
+ {
+ if ((E->trace.rtrac[j][kk]=(double *)malloc(E->trace.itracsize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(initialize tracer arrays)-no memory 1b.%d\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
+ {
+ if ((E->trace.etrac[j][kk]=(double *)malloc(E->trace.itracsize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(initialize tracer arrays)-no memory 1c.%d\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ fprintf(E->trace.fpt,"Physical size of tracer arrays (itracsize): %d\n",
+ E->trace.itracsize[j]);
+ fflush(E->trace.fpt);
+
+return;
+}
+
+
+
+/************ FIND TRACERS *************************************/
+/* */
+/* This function finds tracer elements and moves tracers to */
+/* other processor domains if necessary. */
+/* Array itrac is filled with elemental values. */
+
+void find_tracers(E)
+ struct All_variables *E;
+
+{
+
+int iel;
+int kk;
+int j;
+int it;
+int iprevious_element;
+int num_tracers;
+
+double theta,phi,rad;
+double x,y,z;
+double time_stat1;
+double time_stat2;
+
+int iget_element();
+void put_away_later();
+void eject_tracer();
+void reduce_tracer_arrays();
+void lost_souls();
+void sphere_to_cart();
+
+static int been_here=0;
+
+time_stat1=CPU_time0();
+
+
+if (been_here==0)
+{
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->trace.itrac[j][0];kk++)
+ {
+ E->trace.itrac[j][kk]=-99;
+ }
+ }
+been_here++;
+}
+
+
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+
+/* initialize arrays and statistical counters */
+
+ E->trace.ilater[j]=0;
+
+ E->trace.istat1=0;
+ for (kk=0;kk<=4;kk++)
+ {
+ E->trace.istat_ichoice[j][kk]=0;
+ }
+
+/* important to index by it, not kk */
+
+ it=0;
+ num_tracers=E->trace.itrac[j][0];
+
+ for (kk=1;kk<=num_tracers;kk++)
+ {
+
+ it++;
+
+ theta=E->trace.rtrac[j][0][it];
+ phi=E->trace.rtrac[j][1][it];
+ rad=E->trace.rtrac[j][2][it];
+ x=E->trace.rtrac[j][3][it];
+ y=E->trace.rtrac[j][4][it];
+ z=E->trace.rtrac[j][5][it];
+
+ iprevious_element=E->trace.itrac[j][it];
+
+/* AKMA REMOVE */
+/*
+fprintf(E->trace.fpt,"BB. kk %d %d %d %f %f %f %f %f %f\n",kk,j,iprevious_element,x,y,z,theta,phi,rad);
+fflush(E->trace.fpt);
+*/
+ iel=iget_element(E,j,iprevious_element,x,y,z,theta,phi,rad);
+
+ E->trace.itrac[j][it]=iel;
+
+ if (iel<0)
+ {
+ put_away_later(E,j,it);
+ eject_tracer(E,j,it);
+ it--;
+ }
+
+ } /* end tracers */
+
+} /* end j */
+
+/* Now take care of tracers that exited cap */
+
+/* REMOVE */
+/*
+parallel_process_termination();
+*/
+
+ lost_souls(E);
+
+/* Free later arrays */
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+ if (E->trace.ilater[j]>0)
+ {
+ for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
+ {
+ free(E->trace.rlater[j][kk]);
+ }
+ }
+} /* end j */
+
+
+/* Adjust Array Sizes */
+
+ reduce_tracer_arrays(E);
+
+time_stat2=CPU_time0();
+
+fprintf(E->trace.fpt,"AA: time for find tracers: %f\n", time_stat2-time_stat1);
+
+return;
+}
+
+/************** 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 is receiving from cap n */
+
+
+void lost_souls(E)
+ struct All_variables *E;
+{
+
+
+int ithiscap;
+int ithatcap=1;
+int isend[13][13];
+int ireceive[13][13];
+int isize[13];
+int kk,pp,j;
+int mm;
+int istat_sent[13];
+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 ihorizontal_neighbor;
+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];
+
+int iget_element();
+int icheck_that_processor_shell();
+int icheck_cap();
+void expand_tracer_arrays();
+
+int number_of_caps=12;
+int lev=E->mesh.levmax;
+
+static int been_here=0;
+
+MPI_Status status[200];
+MPI_Request request[200];
+MPI_Status status1;
+MPI_Status status2;
+static int itag=1;
+
+
+/* Note, if for some reason, the number of neighbors exceeds */
+/* 50, which is unlikely, the MPI arrays must be increased. */
+
+
+/* sync all processors */
+
+parallel_process_sync();
+
+/* initialize things */
+
+ if (been_here==0)
+ {
+
+ if (E->parallel.TNUM_PASSz[lev]>2)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-more than 2 zpass\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ if (E->parallel.nprocz>2)
+ {
+ fprintf(stderr,"WARNING(lost souls)-have not tested nprocz>2\n");
+ fprintf(stderr,"WARNING(lost souls)-have not tested nprocz>2\n");
+ fprintf(stderr,"WARNING(lost souls)-have not tested nprocz>2\n");
+ fprintf(E->trace.fpt,"WARNING(lost souls)-have not tested nprocz>2\n");
+ fflush(E->trace.fpt);
+/*
+ if (E->trace.itracer_warnings==1) exit(10);
+*/
+ }
+
+ been_here++;
+ }
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ E->trace.istat_isend=E->trace.ilater[j];
+ }
+
+
+/* initialize isend and ireceive */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ isize[j]=E->trace.ilater[j]*E->trace.number_of_tracer_quantities;
+
+ for (kk=1;kk<=number_of_caps;kk++) isend[j][kk]=0;
+ for (kk=1;kk<=number_of_caps;kk++) ireceive[j][kk]=0;
+ }
+
+/* Allocate Maximum Memory to Send Arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ ithiscap=E->sphere.capid[j];
+
+ itemp_size=max(isize[j],1);
+
+ if ((send[j][ithiscap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+
+ if ((send[j][ithatcap]=(double *)malloc(itemp_size*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+
+/* For testing, remove this */
+/*
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ ithiscap=E->sphere.capid[j];
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+ fprintf(E->trace.fpt,"cap: %d proc %d TNUM: %d ithatcap: %d idd: %d proc_that: %d (%d)\n",
+ ithiscap,E->parallel.me,kk,ithatcap,idd[ithiscap][ithatcap],
+ E->parallel.PROCESSOR[lev][j].pass[kk],E->parallel.mst[ithiscap][ithatcap]);
+
+ }
+ }
+*/
+
+
+/* Pre communication */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+/* transfer tracers from rlater to send */
+
+
+ ithiscap=E->sphere.capid[j];
+
+ numtracers=E->trace.ilater[j];
+
+ for (kk=1;kk<=numtracers;kk++)
+ {
+
+ istat_sent[j]++;
+
+ rad=E->trace.rlater[j][2][kk];
+ x=E->trace.rlater[j][3][kk];
+ y=E->trace.rlater[j][4][kk];
+ z=E->trace.rlater[j][5][kk];
+
+/* first check same cap if nprocz>1 */
+
+ if (E->parallel.nprocz>1)
+ {
+
+ ithatcap=ithiscap;
+ icheck=icheck_cap(E,ithatcap,x,y,z,rad);
+ if (icheck==1) goto foundit;
+
+ }
+
+/* check neighboring caps */
+
+ for (pp=1;pp<=E->parallel.TNUM_PASS[lev][j];pp++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[pp];
+ icheck=icheck_cap(E,ithatcap,x,y,z,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",x,y,z,rad);
+ icheck=icheck_cap(E,ithiscap,x,y,z,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)*E->trace.number_of_tracer_quantities;
+
+ for (pp=0;pp<=(E->trace.number_of_tracer_quantities-1);pp++)
+ {
+ ipos=isend_position+pp;
+ send[j][ithatcap][ipos]=E->trace.rlater[j][pp][kk];
+ }
+
+ } /* end kk, assigning tracers */
+
+ } /* end j */
+
+
+/* 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 */
+/* Also, tags are sent which ensure send and recieve are coordinated. */
+/* E->parallel.mst[i][j]=E->parallel.mst[j][i] */
+
+ idb=0;
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ ithiscap=E->sphere.capid[j];
+
+/* if tracer is in same cap (nprocz>1) */
+
+ if (E->parallel.nprocz>1)
+ {
+ ireceive[j][ithiscap]=isend[j][ithiscap];
+ }
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+
+/* if neighbor cap is in another processor, send information via MPI */
+
+ idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ if (idestination_proc!=E->parallel.me)
+ {
+
+ idb++;
+ MPI_Isend(&isend[j][ithatcap],1,MPI_INT,idestination_proc,
+ E->parallel.mst[ithiscap][ithatcap],MPI_COMM_WORLD,
+ &request[idb-1]);
+
+ } /* end if */
+
+/* if neighbor cap is in the same processor (i.e. less than 12 procs) */
+
+ else if (idestination_proc==E->parallel.me)
+ {
+ for (pp=1;pp<=E->sphere.caps_per_proc;pp++)
+ {
+ if (E->sphere.capid[pp]==ithatcap)
+ {
+ ireceive[pp][ithiscap]=isend[j][ithatcap];
+ goto found_other_cap;
+ }
+ }
+
+ } /* end else if */
+ ;
+found_other_cap:
+ ;
+ } /* end kk, number of neighbors */
+
+ } /* end j, caps per proc */
+
+/* Receive tracer count info */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ ithiscap=E->sphere.capid[j];
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+
+/* if neighbor cap is in another processor, receive information via MPI */
+
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ if (idestination_proc!=E->parallel.me)
+ {
+
+ idb++;
+ MPI_Irecv(&ireceive[j][ithatcap],1,MPI_INT,isource_proc,
+ E->parallel.mst[ithiscap][ithatcap],MPI_COMM_WORLD,
+ &request[idb-1]);
+
+ } /* end if */
+
+ } /* end kk, number of neighbors */
+ } /* end j, caps per proc */
+
+/* Wait for non-blocking calls to complete */
+
+ MPI_Waitall(idb,request,status);
+
+/* Testing, should remove */
+/*
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ ithiscap=E->sphere.capid[j];
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+ fprintf(E->trace.fpt,"cap: %d send %d to cap %d\n",ithiscap,isend[j][ithatcap],ithatcap);
+ fprintf(E->trace.fpt,"cap: %d rec %d from cap %d\n",ithiscap,ireceive[j][ithatcap],ithatcap);
+ }
+ }
+*/
+
+/* Allocate memory in receive arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (ithatcap=1;ithatcap<=number_of_caps;ithatcap++)
+ {
+ isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
+
+ itemp_size=max(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);
+ }
+ }
+ } /* end j */
+
+/* Now, send the tracers to proper caps */
+
+ idb=0;
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ ithiscap=E->sphere.capid[j];
+
+/* same cap */
+
+ if (E->parallel.nprocz>1)
+ {
+
+ ithatcap=ithiscap;
+ isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+ for (mm=0;mm<=(isize[j]-1);mm++)
+ {
+ receive[j][ithatcap][mm]=send[j][ithatcap][mm];
+ }
+
+ }
+
+/* neighbor caps */
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+
+/* if neighbor cap is in another processor, send information via MPI */
+
+ idestination_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ isize[j]=isend[j][ithatcap]*E->trace.number_of_tracer_quantities;
+
+ if (idestination_proc!=E->parallel.me)
+ {
+
+ idb++;
+
+
+ MPI_Isend(send[j][ithatcap],isize[j],MPI_DOUBLE,idestination_proc,
+ E->parallel.mst[ithiscap][ithatcap],MPI_COMM_WORLD,
+ &request[idb-1]);
+
+ } /* end if */
+
+/* if neighbor cap is in the same processor (i.e. less than 12 procs) */
+
+ else if (idestination_proc==E->parallel.me)
+ {
+ for (pp=1;pp<=E->sphere.caps_per_proc;pp++)
+ {
+ if (E->sphere.capid[pp]==ithatcap)
+ {
+
+ for (mm=1;mm<=isize[j];mm++)
+ {
+ receive[pp][ithiscap][mm]=send[j][ithatcap][mm];
+ }
+ goto found_other_cap2;
+ }
+ }
+
+ } /* end else if */
+
+found_other_cap2:
+ ;
+ } /* end kk, number of neighbors */
+
+ } /* end j, caps per proc */
+
+
+/* Receive tracers */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ ithiscap=E->sphere.capid[j];
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+
+/* if neighbor cap is in another processor, receive information via MPI */
+
+ isource_proc=E->parallel.PROCESSOR[lev][j].pass[kk];
+
+ if (idestination_proc!=E->parallel.me)
+ {
+
+ idb++;
+
+ isize[j]=ireceive[j][ithatcap]*E->trace.number_of_tracer_quantities;
+
+ MPI_Irecv(receive[j][ithatcap],isize[j],MPI_DOUBLE,isource_proc,
+ E->parallel.mst[ithiscap][ithatcap],MPI_COMM_WORLD,
+ &request[idb-1]);
+
+ } /* end if */
+
+ } /* end kk, number of neighbors */
+
+ } /* end j, caps per proc */
+
+/* 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) */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ isum[j]=0;
+
+ ithiscap=E->sphere.capid[j];
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+ isum[j]=isum[j]+ireceive[j][ithatcap];
+ }
+ if (E->parallel.nprocz>1) isum[j]=isum[j]+ireceive[j][ithiscap];
+
+ itracers_subject_to_vertical_transport[j]=isum[j];
+ }
+
+/* Allocate Memory for REC array */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(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);
+ }
+ REC[j][0]=0.0;
+ }
+
+/* Put Received tracers in REC */
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ irec[j]=0;
+
+ irec_position=0;
+
+ ithiscap=E->sphere.capid[j];
+
+/* horizontal neighbors */
+
+ for (ihorizontal_neighbor=1;ihorizontal_neighbor<=E->parallel.TNUM_PASS[lev][j];ihorizontal_neighbor++)
+ {
+
+ ithatcap=E->sphere.cap[ithiscap].connection[ihorizontal_neighbor];
+
+ for (pp=1;pp<=ireceive[j][ithatcap];pp++)
+ {
+ irec[j]++;
+ ipos=(pp-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
+ {
+ ipos2=ipos+mm;
+ REC[j][irec_position]=receive[j][ithatcap][ipos2];
+
+ irec_position++;
+
+ } /* end mm (cycling tracer quantities) */
+ } /* end pp (cycling tracers) */
+ } /* end ihorizontal_neighbors (cycling neighbors) */
+
+/* for tracers in the same cap (nprocz>1) */
+
+ if (E->parallel.nprocz>1)
+ {
+ ithatcap=ithiscap;
+ for (pp=1;pp<=ireceive[j][ithatcap];pp++)
+ {
+ irec[j]++;
+ ipos=(pp-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
+ {
+ ipos2=ipos+mm;
+
+ REC[j][irec_position]=receive[j][ithatcap][ipos2];
+
+ irec_position++;
+
+ } /* end mm (cycling tracer quantities) */
+
+ } /* end pp (cycling tracers) */
+
+ } /* endif nproc>1 */
+
+ } /* end j (cycling caps) */
+
+/* Done filling REC */
+
+
+
+/* VERTICAL COMMUNICATION */
+
+/* Note: For generality, I assume that both multiple */
+/* caps per processor as well as multiple processors */
+/* in the radial direction. These are probably */
+/* inconsistent parameter choices for the regular */
+/* CitcomS code. */
+
+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 (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
+ {
+ isize[j]=itracers_subject_to_vertical_transport[j]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(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);
+ }
+ }
+ } /* end j */
+
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+ 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++)
+ {
+
+ it++;
+
+ ireceive_position=((it-1)*E->trace.number_of_tracer_quantities);
+
+ 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_z[j][ivertical_neighbor]++;
+
+ isend_position=(isend_z[j][ivertical_neighbor]-1)*E->trace.number_of_tracer_quantities;
+
+ ilast_receiver_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);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 */
+
+ } /* end j */
+
+
+/* Send arrays are now filled. */
+/* Now send send information to vertical processor neighbor */
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
+
+ MPI_Sendrecv(&isend_z[j][ivertical_neighbor],1,MPI_INT,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag,
+ &ireceive_z[j][ivertical_neighbor],1,MPI_INT,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ itag,MPI_COMM_WORLD,&status1);
+
+/* for testing - remove */
+/*
+fprintf(E->trace.fpt,"PROC: %d IVN: %d (P: %d) SEND: %d REC: %d\n",
+ E->parallel.me,ivertical_neighbor,E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ isend_z[j][ivertical_neighbor],ireceive_z[j][ivertical_neighbor]);
+fflush(E->trace.fpt);
+*/
+
+ } /* end ivertical_neighbor */
+
+} /* end j */
+
+
+/* Allocate memory to receive_z arrays */
+
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->parallel.TNUM_PASSz[lev];kk++)
+ {
+ isize[j]=ireceive_z[j][kk]*E->trace.number_of_tracer_quantities;
+ isize[j]=max(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);
+ }
+ }
+ } /* end j */
+
+/* Send Tracers */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
+ isize_send=isend_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+ isize_receive=ireceive_z[j][ivertical_neighbor]*E->trace.number_of_tracer_quantities;
+
+ MPI_Sendrecv(send_z[j][ivertical_neighbor],isize_send,
+ MPI_DOUBLE,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],itag+1,
+ receive_z[j][ivertical_neighbor],isize_receive,
+ MPI_DOUBLE,
+ E->parallel.PROCESSORz[lev].pass[ivertical_neighbor],
+ itag+1,MPI_COMM_WORLD,&status2);
+
+ }
+ }
+
+/* Put tracers into REC array */
+
+/* First, reallocate memory to REC */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ isum[j]=0;
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
+ isum[j]=isum[j]+ireceive_z[j][ivertical_neighbor];
+ }
+
+ isum[j]=isum[j]+irec[j];
+
+ isize[j]=isum[j]*E->trace.number_of_tracer_quantities;
+
+ 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 (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (ivertical_neighbor=1;ivertical_neighbor<=E->parallel.TNUM_PASSz[lev];ivertical_neighbor++)
+ {
+
+ for (kk=1;kk<=ireceive_z[j][ivertical_neighbor];kk++)
+ {
+ irec[j]++;
+
+ irec_position=(irec[j]-1)*E->trace.number_of_tracer_quantities;
+ ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_tracer_quantities-1);mm++)
+ {
+ REC[j][irec_position+mm]=receive_z[j][ivertical_neighbor][ireceive_position+mm];
+ }
+ }
+
+ }
+ }
+
+/* Free Vertical Arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ 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 (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=irec[j];kk++)
+ {
+ E->trace.itrac[j][0]++;
+
+ if (E->trace.itrac[j][0]>(E->trace.itracsize[j]-5)) expand_tracer_arrays(E,j);
+
+ ireceive_position=(kk-1)*E->trace.number_of_tracer_quantities;
+
+ for (mm=0;mm<=(E->trace.number_of_advection_quantities-1);mm++)
+ {
+ ipos=ireceive_position+mm;
+
+ E->trace.rtrac[j][mm][E->trace.itrac[j][0]]=REC[j][ipos];
+ }
+ for (mm=0;mm<=(E->trace.number_of_extra_quantities-1);mm++)
+ {
+ ipos=ireceive_position+E->trace.number_of_advection_quantities+mm;
+
+ E->trace.etrac[j][mm][E->trace.itrac[j][0]]=REC[j][ipos];
+ }
+
+ theta=E->trace.rtrac[j][0][E->trace.itrac[j][0]];
+ phi=E->trace.rtrac[j][1][E->trace.itrac[j][0]];
+ rad=E->trace.rtrac[j][2][E->trace.itrac[j][0]];
+ x=E->trace.rtrac[j][3][E->trace.itrac[j][0]];
+ y=E->trace.rtrac[j][4][E->trace.itrac[j][0]];
+ z=E->trace.rtrac[j][5][E->trace.itrac[j][0]];
+
+
+ iel=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);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ E->trace.itrac[j][E->trace.itrac[j][0]]=iel;
+
+ }
+ }
+
+
+/* Free Arrays */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ ithiscap=E->sphere.capid[j];
+
+ free(REC[j]);
+
+ free(send[j][ithiscap]);
+
+ for (kk=1;kk<=E->parallel.TNUM_PASS[lev][j];kk++)
+ {
+ ithatcap=E->sphere.cap[ithiscap].connection[kk];
+
+ free(send[j][ithatcap]);
+
+ }
+
+ for (ithatcap=1;ithatcap<=number_of_caps;ithatcap++)
+ {
+ free(receive[j][ithatcap]);
+ }
+
+ }
+
+return;
+}
+
+
+/****** REDUCE TRACER ARRAYS *****************************************/
+
+void reduce_tracer_arrays(E)
+ struct All_variables *E;
+
+{
+
+int inewsize;
+int kk;
+int iempty_space;
+int j;
+
+int icushion=100;
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+
+/* if physical size is double tracer size, reduce it */
+
+ iempty_space=(E->trace.itracsize[j]-E->trace.itrac[j][0]);
+
+ if (iempty_space>(E->trace.itrac[j][0]+icushion))
+ {
+
+
+ inewsize=E->trace.itrac[j][0]+E->trace.itrac[j][0]/4+icushion;
+
+ if (inewsize<1)
+ {
+ fprintf(E->trace.fpt,"Error(reduce tracer arrays)-something up (hdf3)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+ if ((E->trace.itrac[j]=(int *)realloc(E->trace.itrac[j],inewsize*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(reduce tracer arrays )-no memory (itracer)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+ for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
+ {
+ if ((E->trace.rtrac[j][kk]=(double *)realloc(E->trace.rtrac[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
+ {
+ if ((E->trace.etrac[j][kk]=(double *)realloc(E->trace.etrac[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(reduce tracer arrays )-no memory 783 (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+
+ fprintf(E->trace.fpt,"Reducing physical memory of itrac, rtrac, and etrac to %d from %d\n",
+ E->trace.itracsize[j],inewsize);
+
+ E->trace.itracsize[j]=inewsize;
+
+ } /* end if */
+
+} /* end j */
+
+return;
+}
+
+/********** PUT AWAY LATER ************************************/
+/* */
+/* rlater has a similar structure to rtrac */
+/* ilatersize is the physical memory and */
+/* ilater is the number of tracers */
+
+
+void put_away_later(E,j,it)
+ struct All_variables *E;
+ int it,j;
+
+{
+
+
+int kk;
+
+
+void expand_later_array();
+
+
+/* The first tracer in initiates memory allocation. */
+/* Memory is freed after parallel communications */
+
+ if (E->trace.ilater[j]==0)
+ {
+
+ E->trace.ilatersize[j]=E->trace.itracsize[j]/5;
+
+ for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
+ {
+ if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(put_away_later)-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ } /* end first particle initiating memory allocation */
+
+
+/* Put tracer in later array */
+
+ E->trace.ilater[j]++;
+
+ if (E->trace.ilater[j] >= (E->trace.ilatersize[j]-5)) expand_later_array(E,j);
+
+/* stack advection and extra quantities together (advection first) */
+
+ for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
+ {
+ E->trace.rlater[j][kk][E->trace.ilater[j]]=E->trace.rtrac[j][kk][it];
+ }
+ for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
+ {
+ E->trace.rlater[j][E->trace.number_of_advection_quantities+kk][E->trace.ilater[j]]=E->trace.etrac[j][kk][it];
+ }
+
+return;
+}
+
+/****** EXPAND LATER ARRAY *****************************************/
+
+void expand_later_array(E,j)
+ struct All_variables *E;
+ int j;
+{
+
+int inewsize;
+int kk;
+int icushion;
+
+/* expand rlater by 20% */
+
+ icushion=100;
+
+ inewsize=E->trace.ilatersize[j]+E->trace.ilatersize[j]/5+icushion;
+
+ for (kk=0;kk<=((E->trace.number_of_tracer_quantities)-1);kk++)
+ {
+ if ((E->trace.rlater[j][kk]=(double *)realloc(E->trace.rlater[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(expand later array )-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+
+ fprintf(E->trace.fpt,"Expanding physical memory of rlater to %d from %d\n",
+ inewsize,E->trace.ilatersize[j]);
+
+ E->trace.ilatersize[j]=inewsize;
+
+return;
+}
+
+
+/***** EJECT TRACER ************************************************/
+
+
+void eject_tracer(E,j,it)
+ struct All_variables *E;
+ int it,j;
+
+{
+
+int ilast_tracer;
+int kk;
+
+
+ ilast_tracer=E->trace.itrac[j][0];
+
+/* put last tracer in ejected tracer position */
+
+ E->trace.itrac[j][it]=E->trace.itrac[j][ilast_tracer];
+
+ for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
+ {
+ E->trace.rtrac[j][kk][it]=E->trace.rtrac[j][kk][ilast_tracer];
+ }
+ for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
+ {
+ E->trace.etrac[j][kk][it]=E->trace.etrac[j][kk][ilast_tracer];
+ }
+
+
+ E->trace.itrac[j][0]--;
+
+return;
+}
+
+
+
+/*********** MAKE REGULAR GRID ********************************/
+/* */
+/* This function generates the finer regular grid which is */
+/* mapped to real elements */
+
+void make_regular_grid(E)
+ 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;
+
+int icheck_all_columns();
+int icheck_element_column();
+int icheck_column_neighbors();
+void sphere_to_cart();
+void fix_theta();
+void fix_phi();
+
+static double two_pi=2.0*M_PI;
+
+ elz=E->lmesh.elz;
+
+ nelsurf=E->lmesh.elx*E->lmesh.ely;
+
+/* 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");
+ fflush(stderr);
+
+
+/* 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=max(thetamax,theta);
+ thetamin=min(thetamin,theta);
+
+ }
+
+/* expand range slightly (should take care of poles) */
+
+ thetamax=thetamax+expansion;
+ thetamax=min(thetamax,M_PI);
+
+ thetamin=thetamin-expansion;
+ thetamin=max(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);
+ fflush(stderr);
+ if (E->trace.itracer_warnings==1) 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 (-99 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 =-99 */
+
+ for (kk=1;kk<=numregnodes;kk++)
+ {
+ E->trace.regnodetoel[j][kk]=-99;
+ }
+
+/* Begin Mapping (only need to use surface elements) */
+
+ parallel_process_sync();
+ if (E->parallel.me==0) fprintf(stderr,"Beginning Mapping\n");
+ fflush(stderr);
+
+/* 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=min(theta_min,theta);
+ theta_max=max(theta_max,theta);
+ phi_min=min(phi_min,phi);
+ phi_max=max(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=min(M_PI,theta_max);
+ theta_min=max(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_phi(&phi_max);
+ fix_phi(&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]=-99;
+
+/* 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);
+
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+
+
+ 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,x,y,z,rad);
+ if (ival>0)
+ {
+ E->trace.regnodetoel[j][kk]=ilast_el;
+ goto foundit;
+ }
+
+/* check neighbors */
+
+ ival=icheck_column_neighbors(E,j,ilast_el,x,y,z,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,x,y,z,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]!=-99)
+ {
+ 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");
+ fflush(stderr);
+
+/* 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 = -99 (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 (-99). */
+/* 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();
+
+ if (E->parallel.me==0) fprintf(stderr,"Mapping completed (%f seconds)\n",CPU_time0()-start_time);
+ fflush(stderr);
+
+/* 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=1 is optimal)\n");
+ fprintf(E->trace.fpt," (if ichoice=0, no elements are touched by regular element)\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 */
+
+
+return;
+}
+
+
+
+/********* ICHECK COLUMN NEIGHBORS ***************************/
+/* */
+/* This function check whether a point is in a neighboring */
+/* column. Neighbor surface element number is returned */
+
+int icheck_column_neighbors(E,j,nel,x,y,z,rad)
+ struct All_variables *E;
+ int j,nel;
+ double x,y,z,rad;
+{
+
+int ival;
+int neighbor[25];
+int elx,ely,elz;
+int elxz;
+int kk;
+
+int icheck_element_column();
+
+/*
+static int number_of_neighbors=24;
+*/
+
+/* maybe faster to only check inner ring */
+
+static 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],x,y,z,rad);
+ if (ival>0)
+ {
+ return neighbor[kk];
+ }
+ }
+ }
+
+return -99;
+}
+
+/********** 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 -99 if can't be found. */
+
+int icheck_all_columns(E,j,x,y,z,rad)
+ struct All_variables *E;
+ int j;
+ double x,y,z,rad;
+{
+
+int icheck;
+int nel;
+int icheck_element_column();
+
+int elz=E->lmesh.elz;
+int numel=E->lmesh.nel;
+
+ for (nel=elz;nel<=numel;nel=nel+elz)
+ {
+ icheck=icheck_element_column(E,j,nel,x,y,z,rad);
+ if (icheck==1)
+ {
+ return nel;
+ }
+ }
+
+
+return -99;
+}
+
+
+/************ INITIALIZE OLD COMPOSITION ************************/
+void initialize_old_composition(E)
+ struct All_variables *E;
+{
+
+char output_file[200];
+char input_s[1000];
+
+int ibottom_node;
+int kk;
+int istep;
+int idum;
+int j;
+
+double zbottom;
+double time;
+
+FILE *fp;
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+ if ((E->trace.oldel[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(fill old composition)-no memory 324c\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+}
+
+
+ if ((E->trace.itracer_restart==0)||(E->trace.itracer_restart==-1))
+ {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+
+ ibottom_node=E->ien[j][kk].node[1];
+ zbottom=E->sx[j][3][ibottom_node];
+
+ if (zbottom<E->trace.z_interface) E->trace.oldel[j][kk]=1.0;
+ if (zbottom>=E->trace.z_interface) E->trace.oldel[j][kk]=0.0;
+
+ } /* end kk */
+ } /* end j */
+ }
+
+
+/* Else read from file */
+
+
+ else if (E->trace.itracer_restart==1)
+ {
+
+/* first look for backing file */
+
+
+ sprintf(output_file,"%s.OLDCOMP.%d.%d",E->control.old_P_file,E->parallel.me,E->control.restart);
+ if ( (fp=fopen(output_file,"r"))==NULL)
+ {
+ sprintf(output_file,"%s.OLDCOMP.%d",E->control.old_P_file,E->parallel.me);
+
+ if ( (fp=fopen(output_file,"r"))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKMerror(Initialize Old Composition)-OLDCOMP NOT EXIST\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ fgets(input_s,200,fp);
+ sscanf(input_s,"%d %d %lf",&idum,&istep,&time);
+
+ if (istep!=E->control.restart)
+ {
+ fprintf(E->trace.fpt,"Error(initialize old comp) %d %d\n",
+ istep,E->monitor.solution_cycles);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ fgets(input_s,200,fp);
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ fgets(input_s,200,fp);
+ sscanf(input_s,"%lf",&E->trace.oldel[j][kk]);
+ }
+ }
+
+ fclose(fp);
+
+ } /* endif */
+
+
+
+return;
+}
+
+/************ FILL COMPOSITION ************************/
+void fill_composition(E)
+ struct All_variables *E;
+{
+
+void compute_elemental_composition_ratio_method();
+void map_composition_to_nodes();
+
+/* Currently, only the ratio method works here. */
+/* Will have to come back here to include the absolute method. */
+
+/* ratio method */
+
+ if (E->trace.ibuoy_type==1)
+ {
+ compute_elemental_composition_ratio_method(E);
+ }
+
+/* absolute method */
+
+ if (E->trace.ibuoy_type!=1)
+ {
+ fprintf(E->trace.fpt,"Error(compute...)-only ratio method now\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+/* Map elemental composition to nodal points */
+
+ map_composition_to_nodes(E);
+
+return;
+}
+
+/*********** COMPUTE ELEMENTAL COMPOSITION RATIO METHOD ***/
+/* */
+/* This function computes the composition per element. */
+/* This function computes the composition per element. */
+/* Integer array ieltrac stores tracers per element. */
+/* Double array celtrac stores the sum of tracer composition */
+
+void compute_elemental_composition_ratio_method(E)
+ struct All_variables *E;
+{
+
+int kk;
+int numtracers;
+int nelem;
+int j;
+int iempty=0;
+
+double comp;
+
+static int been_here=0;
+
+if (been_here==0)
+{
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ if ((E->trace.ieltrac[j]=(int *)malloc((E->lmesh.nel+1)*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 5u83a\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((E->trace.celtrac[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 58hy8\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ if ((E->trace.comp_el[j]=(double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(compute_elemental_composition)-no memory 8989y\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ been_here++;
+
+}
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+/* first zero arrays */
+
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ E->trace.ieltrac[j][kk]=0;
+ E->trace.celtrac[j][kk]=0.0;
+ }
+
+ numtracers=E->trace.itrac[j][0];
+
+/* Fill ieltrac and celtrac */
+
+
+ for (kk=1;kk<=numtracers;kk++)
+ {
+
+ nelem=E->trace.itrac[j][kk];
+ E->trace.ieltrac[j][nelem]++;
+
+ comp=E->trace.etrac[j][0][kk];
+
+ if (comp>1.0000001)
+ {
+ fprintf(E->trace.fpt,"ERROR(compute elemental)-not ready for comp>1 yet (%f)(tr. %d) \n",comp,kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ E->trace.celtrac[j][nelem]=E->trace.celtrac[j][nelem]+comp;
+
+ }
+
+/* Check for empty entries and compute ratio. */
+/* If no tracers are in an element, use previous composition */
+
+ iempty=0;
+
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+
+ if (E->trace.ieltrac[j][kk]==0)
+ {
+ iempty++;
+ E->trace.comp_el[j][kk]=E->trace.oldel[j][kk];
+ }
+ else if (E->trace.ieltrac[j][kk]>0)
+ {
+ E->trace.comp_el[j][kk]=E->trace.celtrac[j][kk]/(1.0*E->trace.ieltrac[j][kk]);
+ }
+
+ if (E->trace.comp_el[j][kk]>(1.000001) || E->trace.comp_el[j][kk]<(-0.000001))
+ {
+ fprintf(E->trace.fpt,"ERROR(compute elemental)-noway (3u5hd)\n");
+ fprintf(E->trace.fpt,"COMPEL: %f (%d)(%d)\n",E->trace.comp_el[j][kk],kk,E->trace.ieltrac[j][kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ if (iempty>0)
+ {
+
+/*
+ fprintf(E->trace.fpt,"%d empty elements filled with old values (%f percent)\n",iempty, (100.0*iempty/E->lmesh.nel));
+ fflush(E->trace.fpt);
+*/
+
+ if ((1.0*iempty/E->lmesh.nel)>0.80)
+ {
+ fprintf(E->trace.fpt,"WARNING(compute_elemental...)-number of tracers is REALLY LOW\n");
+ fflush(E->trace.fpt);
+ if (E->trace.itracer_warnings==1) exit(10);
+ }
+ }
+
+/* Fill oldel */
+
+
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ E->trace.oldel[j][kk]=E->trace.comp_el[j][kk];
+ }
+
+} /* end j */
+
+ E->trace.istat_iempty=E->trace.istat_iempty+iempty;
+
+return;
+}
+
+/********** MAP COMPOSITION TO NODES ****************/
+/* */
+
+
+void map_composition_to_nodes(E)
+ struct All_variables *E;
+{
+
+int kk;
+int nelem, nodenum;
+int j;
+
+void exchange_node_d();
+
+static int been_here=0;
+
+if (been_here==0)
+{
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ if ((E->trace.comp_node[j]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"AKM(map_compostion_to_nodes)-no memory 983rk\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+been_here++;
+}
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+/* first, initialize node array */
+
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ {
+ E->trace.comp_node[j][kk]=0.0;
+ }
+
+/* Loop through all elements */
+
+ for (nelem=1;nelem<=E->lmesh.nel;nelem++)
+ {
+
+/* for each element, loop through element nodes */
+
+/* weight composition */
+
+ for (nodenum=1;nodenum<=8;nodenum++)
+ {
+
+ E->trace.comp_node[j][E->ien[j][nelem].node[nodenum]] +=
+ E->trace.comp_el[j][nelem]*
+ E->TWW[E->mesh.levmax][j][nelem].node[nodenum];
+
+ }
+
+ } /* end nelem */
+} /* end j */
+
+/* akm modified exchange node routine for doubles */
+
+ exchange_node_d(E,E->trace.comp_node,E->mesh.levmax);
+
+/* Divide by nodal volume */
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ {
+ E->trace.comp_node[j][kk] *= E->MASS[E->mesh.levmax][j][kk];
+ }
+
+/* testing */
+/*
+ for (kk=1;kk<=E->lmesh.nel;kk++)
+ {
+ fprintf(E->trace.fpt,"%d %f\n",kk,E->trace.comp_el[j][kk]);
+ }
+
+ for (kk=1;kk<=E->lmesh.nno;kk++)
+ {
+ fprintf(E->trace.fpt,"%d %f %f\n",kk,E->sx[j][3][kk],E->trace.comp_node[j][kk]);
+ }
+*/
+
+
+} /* end j */
+
+return;
+}
+
+/**** WRITE TRACE INSTRUCTIONS ***************/
+void write_trace_instructions(E)
+ struct All_variables *E;
+{
+
+
+
+
+ fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
+ fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
+
+ if (E->trace.itracer_type==0)
+ {
+ fprintf(E->trace.fpt,"Passive Tracers\n");
+ }
+
+ if (E->trace.itracer_type==1)
+ {
+ fprintf(E->trace.fpt,"Active Tracers\n");
+ if (E->trace.ibuoy_type==1) fprintf(E->trace.fpt,"Ratio Method\n");
+ if (E->trace.ibuoy_type==0) fprintf(E->trace.fpt,"Absolute Method\n");
+ fprintf(E->trace.fpt,"Buoyancy Ratio: %f\n", E->trace.buoyancy_ratio);
+
+ if (E->trace.ireset_initial_composition==0)
+ {
+ fprintf(E->trace.fpt,"Using old initial composition from tracer files\n");
+ }
+ else
+ {
+ fprintf(E->trace.fpt,"Resetting initial composition\n");
+ }
+
+ if (E->trace.itrace_filter==0)
+ {
+ fprintf(E->trace.fpt,"Filter is OFF\n");
+ }
+ else if (E->trace.itrace_filter==1)
+ {
+ fprintf(stderr,"Filter is ON! - Type 1\n");
+ fprintf(E->trace.fpt,"Filter is ON! - Type 1\n");
+ fprintf(E->trace.fpt," Compositions greater than %f are increased\n",E->trace.trace_filter_uplimit);
+ fprintf(E->trace.fpt," Compositions less than %f are decreased\n",E->trace.trace_filter_downlimit);
+ }
+ else if (E->trace.itrace_filter==2)
+ {
+ fprintf(stderr,"Filter is ON! - Type 2\n");
+ fprintf(E->trace.fpt,"Filter is ON! - Type 2\n");
+ }
+ else
+ {
+ fprintf(stderr,"Invalid Filter type %d\n",E->trace.itrace_filter);
+ fprintf(E->trace.fpt,"Invalid Filter type %d\n",E->trace.itrace_filter);
+ }
+ fflush(E->trace.fpt);
+ fflush(stderr);
+
+ }
+
+ if (E->trace.itracer_restart==0)
+ {
+ fprintf(E->trace.fpt,"Generating New Tracer Array\n");
+ fprintf(E->trace.fpt,"Tracers per element: %d\n",E->trace.itperel);
+ if (E->trace.itracer_type==1)
+ {
+ fprintf(E->trace.fpt,"Interface Height: %f\n",E->trace.z_interface);
+ }
+ }
+ if (E->trace.itracer_restart==-1)
+ {
+ fprintf(E->trace.fpt,"Reading tracer file %s\n",E->trace.tracer_file);
+ if (E->trace.icartesian_or_spherical_input==0)
+ {
+ fprintf(E->trace.fpt,"Coordinates are read as Cartesian\n");
+ fprintf(E->trace.fpt," x,y,z format\n");
+ }
+ else if (E->trace.icartesian_or_spherical_input==1)
+ {
+ fprintf(E->trace.fpt,"Coordinates are read as Spherical\n");
+ fprintf(E->trace.fpt," theta,phi,rad format\n");
+ }
+ else
+ {
+ fprintf(E->trace.fpt,"ERROR-Cartesian or Spherical input ?\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ }
+ if (E->trace.itracer_restart==1)
+ {
+ fprintf(E->trace.fpt,"Restarting Tracers\n");
+ }
+
+ if (E->trace.itracer_advection_scheme==1)
+ {
+ fprintf(E->trace.fpt,"\nSimple predictor-corrector method used\n");
+ fprintf(E->trace.fpt,"(Uses only velocity at to) \n");
+ fprintf(E->trace.fpt,"(xf=x0+0.5*dt*(v(x0,y0,z0) + v(xp,yp,zp))\n\n");
+ }
+ else if (E->trace.itracer_advection_scheme==2)
+ {
+ fprintf(E->trace.fpt,"\nPredictor-corrector method used\n");
+ fprintf(E->trace.fpt,"(Uses only velocity at to and to+dt) \n");
+ fprintf(E->trace.fpt,"(xf=x0+0.5*dt*(v(x0,y0,z0,to) + v(xp,yp,zp,to+dt))\n\n");
+ }
+ else
+ {
+ fprintf(E->trace.fpt,"Sorry-Other Advection Schemes are Unavailable %d\n",E->trace.itracer_advection_scheme);
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+ if (E->trace.itracer_interpolation_scheme==1)
+ {
+ fprintf(E->trace.fpt,"\nGreat Circle Projection Interpolation Scheme \n");
+ }
+ else if (E->trace.itracer_interpolation_scheme==2)
+ {
+ fprintf(E->trace.fpt,"\nSimple Averaging Interpolation Scheme \n");
+ fprintf(E->trace.fpt,"\n(Not that great of a scheme!) \n");
+
+ fprintf(E->trace.fpt,"Sorry-Other Interpolation Schemes are Unavailable\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+
+ }
+ else
+ {
+ fprintf(E->trace.fpt,"Sorry-Other Interpolation Schemes are Unavailable\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+ fprintf(E->trace.fpt,"Writing Tracers to File every %d steps\n",
+ E->trace.iwrite_tracers_every);
+
+
+ if (E->trace.icompositional_rheology==0)
+ {
+ fprintf(E->trace.fpt,"Compositional Rheology - OFF\n");
+ }
+ else if (E->trace.icompositional_rheology>0)
+ {
+ fprintf(E->trace.fpt,"Compositional Rheology - ON\n");
+ fprintf(E->trace.fpt,"Compositional Prefactor: %f\n",
+ E->trace.compositional_rheology_prefactor);
+ }
+
+/* regular grid stuff */
+
+ fprintf(E->trace.fpt,"Regular Grid-> deltheta: %f delphi: %f\n",
+ E->trace.deltheta[0],E->trace.delphi[0]);
+
+
+
+
+/* more obscure stuff */
+
+ fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
+ fprintf(E->trace.fpt,"Number of Advection Quantities: %d\n",
+ E->trace.number_of_advection_quantities);
+ fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
+ E->trace.number_of_extra_quantities);
+ fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
+ E->trace.number_of_tracer_quantities);
+
+
+/* analytical test */
+
+ if (E->trace.ianalytical_tracer_test==1)
+ {
+ fprintf(E->trace.fpt,"\n\n ! Analytical Test Being Performed ! \n");
+ fprintf(E->trace.fpt,"(some of the above parameters may not be used or applied\n");
+ fprintf(E->trace.fpt,"Velocity functions given in main code\n");
+ fflush(E->trace.fpt);
+ }
+
+ if (E->trace.itracer_warnings==0)
+ {
+ fprintf(E->trace.fpt,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
+ fprintf(stderr,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ }
+
+return;
+}
+
+/******* VISCOSITY CHECKS *********************/
+void viscosity_checks(E)
+ struct All_variables *E;
+
+{
+
+int kk;
+
+
+ fprintf(E->trace.fpt,"Sorry, compositional rheology not yet available\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+
+ if (E->viscosity.update_allowed==0)
+ {
+ for (kk=1;kk<=5;kk++)
+ {
+ if (E->parallel.me==0)
+ {
+ fprintf(stdout,"I hope you know viscosity update is off!\n");
+ fprintf(stdout," This better be a fixed boundary!\n");
+ }
+
+ fprintf(E->trace.fpt,"\nI hope you know viscosity update is off!\n");
+ fprintf(E->trace.fpt," This better be a fixed boundary!\n\n");
+ fflush(stdout);
+ fflush(E->trace.fpt);
+ }
+ }
+ if (E->viscosity.TDEPV!=0)
+ {
+ fprintf(E->trace.fpt,"Sorry, must fix comp rheology for TDEPV\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+ if (E->viscosity.SDEPV!=0)
+ {
+ fprintf(E->trace.fpt,"Sorry, must fix comp rheology for SDEPV\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+
+return;
+}
+
+/************** SETUP SHARED CAP INFORMATION ***************/
+void setup_shared_cap_information(E)
+ struct All_variables *E;
+{
+
+int node[5];
+int nox,noy,noz;
+int ithiscap;
+int j;
+int kk;
+int nroot;
+int nproc;
+int irootcap;
+
+double x,y,z;
+double theta,phi,rad;
+double diff;
+
+static double eps=0.0001;
+
+void sphere_to_cart();
+
+
+ nox=E->lmesh.nox;
+ noy=E->lmesh.noy;
+ noz=E->lmesh.noz;
+
+ node[1]=nox*noz*(noy-1)+noz;
+ node[2]=noz;
+ node[3]=noz*nox;
+ node[4]=noz*nox*noy;
+
+/* determine coords of cap surface nodes */
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+ ithiscap=E->sphere.capid[j];
+
+ E->trace.XCAP[ithiscap][0]=0.0;
+ E->trace.YCAP[ithiscap][0]=0.0;
+ E->trace.ZCAP[ithiscap][0]=0.0;
+ E->trace.THETA_CAP[ithiscap][0]=0.0;
+ E->trace.PHI_CAP[ithiscap][0]=0.0;
+ E->trace.RAD_CAP[ithiscap][0]=0.0;
+
+ for (kk=1;kk<=4;kk++)
+ {
+
+ theta=E->sx[j][1][node[kk]];
+ phi=E->sx[j][2][node[kk]];
+ rad=E->sphere.ro;
+
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+
+ E->trace.XCAP[ithiscap][kk]=x;
+ E->trace.YCAP[ithiscap][kk]=y;
+ E->trace.ZCAP[ithiscap][kk]=z;
+ E->trace.THETA_CAP[ithiscap][kk]=theta;
+ E->trace.PHI_CAP[ithiscap][kk]=phi;
+ E->trace.RAD_CAP[ithiscap][kk]=rad;
+ E->trace.COS_THETA[ithiscap][kk]=cos(theta);
+ E->trace.SIN_THETA[ithiscap][kk]=sin(theta);
+ E->trace.COS_PHI[ithiscap][kk]=cos(phi);
+ E->trace.SIN_PHI[ithiscap][kk]=sin(phi);
+
+ }
+
+}
+
+
+/* broadcast information */
+
+ parallel_process_sync();
+ nproc=E->parallel.nproc;
+
+ for (nroot=0;nroot<=(nproc-1);nroot++)
+ {
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+
+ irootcap=E->sphere.capid[j];
+
+
+ MPI_Bcast(&irootcap,1,MPI_INT,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.XCAP[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.YCAP[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.ZCAP[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.THETA_CAP[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.PHI_CAP[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.RAD_CAP[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.COS_THETA[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.SIN_THETA[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.COS_PHI[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(E->trace.SIN_PHI[irootcap],5,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ }
+ }
+
+/* Share some processor information */
+
+ parallel_process_sync();
+
+ if (nproc>100)
+ {
+ fprintf(E->trace.fpt,"Error(setup_shared_cap)-should increase arrays TOP AND BOTTOM_RAD\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ E->trace.BOTTOM_RAD[E->parallel.me][j]=E->sx[j][3][1];
+ E->trace.TOP_RAD[E->parallel.me][j]=E->sx[j][3][noz];
+
+ diff=fabs(E->sx[j][3][noz]-E->sphere.ro);
+
+ if (diff<=eps) E->trace.ITOP[E->parallel.me][j]=1;
+ else E->trace.ITOP[E->parallel.me][j]=0;
+
+
+ }
+
+/* broadcast */
+
+ for (nroot=0;nroot<=(nproc-1);nroot++)
+ {
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+
+ MPI_Bcast(&E->trace.BOTTOM_RAD[nroot][j],1,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+ MPI_Bcast(&E->trace.TOP_RAD[nroot][j],1,MPI_DOUBLE,
+ nroot,MPI_COMM_WORLD);
+
+
+ }
+ }
+
+return;
+}
+
+/************** RESTART TRACERS ******************************************/
+/* */
+/* This function restarts tracers written from previous calculation */
+/* and the tracers are read as seperate files for each processor domain. */
+
+void restart_tracers(E)
+ struct All_variables *E;
+{
+
+char output_file[200];
+char input_s[1000];
+
+int j,kk;
+int idum1,istep;
+int numtracers;
+
+double rdum1;
+double theta,phi,rad;
+double comp;
+double x,y,z;
+
+void initialize_tracer_arrays();
+void sphere_to_cart();
+void keep_in_sphere();
+
+FILE *fp1;
+
+
+/* first try old storage file */
+
+ sprintf(output_file,"%s.tracers.%d.%d",E->control.old_P_file,E->parallel.me,E->control.restart);
+
+/* if not, try last backing file */
+
+ if ( (fp1=fopen(output_file,"r"))==NULL)
+ {
+ sprintf(output_file,"%s.tracers.%d",E->control.old_P_file,E->parallel.me);
+ if ( (fp1=fopen(output_file,"r"))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-file not found %s\n",output_file);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ fprintf(stderr,"Restarting Tracers from %s\n",output_file);
+ fflush(stderr);
+
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ fgets(input_s,200,fp1);
+ sscanf(input_s,"%d %d %lf %lf %lf %d",
+ &numtracers,&istep,&rdum1,&E->trace.bulk_composition,&E->trace.initial_bulk_composition,&idum1);
+
+
+
+/* some error control */
+
+ if (istep!=E->control.restart)
+ {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)- step? %d %d\n",istep,E->control.restart);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+/* allocate memory for tracer arrays */
+
+ initialize_tracer_arrays(E,j,numtracers);
+ E->trace.itrac[j][0]=numtracers;
+
+ for (kk=1;kk<=numtracers;kk++)
+ {
+ comp=0.0;
+ fgets(input_s,200,fp1);
+ if (E->trace.itracer_type==0)
+ {
+ sscanf(input_s,"%lf %lf %lf\n",&theta,&phi,&rad);
+ }
+ else if (E->trace.itracer_type==1)
+ {
+ sscanf(input_s,"%lf %lf %lf %lf\n",&theta,&phi,&rad,&comp);
+ }
+ else
+ {
+ fprintf(E->trace.fpt,"ERROR(restart tracers)-huh?\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ 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 */
+
+ keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+
+ E->trace.rtrac[j][0][kk]=theta;
+ E->trace.rtrac[j][1][kk]=phi;
+ E->trace.rtrac[j][2][kk]=rad;
+ E->trace.rtrac[j][3][kk]=x;
+ E->trace.rtrac[j][4][kk]=y;
+ E->trace.rtrac[j][5][kk]=z;
+ if (E->trace.itracer_type==1) E->trace.etrac[j][0][kk]=comp;
+
+ }
+
+ fprintf(E->trace.fpt,"Read %d tracers from file %s\n",numtracers,output_file);
+ fflush(E->trace.fpt);
+
+ }
+
+
+return;
+}
+
+/************** MAKE TRACER ARRAY ********************************/
+/* Here, each cap will generate tracers somewhere */
+/* in the sphere - check if its in this cap - then check radial */
+
+void make_tracer_array(E)
+ struct All_variables *E;
+{
+
+int kk;
+int tracers_cap;
+int j;
+int ithiscap;
+int ival;
+int number_of_tries=0;
+int max_tries;
+
+double x,y,z;
+double theta,phi,rad;
+double dmin,dmax;
+double random1,random2,random3;
+double processor_fraction;
+
+void cart_to_sphere();
+void keep_in_sphere();
+void initialize_tracer_arrays();
+
+int icheck_cap();
+
+if (E->parallel.me==0) fprintf(stderr,"Making Tracer Array\n");
+fflush(stderr);
+
+
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+
+ ithiscap=E->sphere.capid[j];
+
+
+ processor_fraction=( ( pow(E->sx[j][3][E->lmesh.noz],3.0)-pow(E->sx[j][3][1],3.0))/
+ (pow(E->sphere.ro,3.0)-pow(E->sphere.ri,3.0)));
+
+ tracers_cap=((E->lmesh.nel*E->parallel.nprocz)*E->trace.itperel
+ *processor_fraction);
+/*
+fprintf(stderr,"AA: proc frac: %f (%d) %d %d %f %f\n",processor_fraction,tracers_cap,E->lmesh.nel,E->parallel.nprocz, E->sx[j][3][E->lmesh.noz],E->sx[j][3][1]);
+*/
+
+ fprintf(E->trace.fpt,"\nGenerating %d Tracers\n",tracers_cap);
+
+ initialize_tracer_arrays(E,j,tracers_cap);
+
+
+/* Tracers are placed randomly in cap */
+/* (intentionally using rand() instead of srand() )*/
+
+ dmin=-1.0*E->sphere.ro;
+ dmax=E->sphere.ro;
+
+ while (E->trace.itrac[j][0]<tracers_cap)
+ {
+
+ number_of_tries++;
+ max_tries=500*tracers_cap;
+
+ if (number_of_tries>max_tries)
+ {
+ fprintf(E->trace.fpt,"Error(make_tracer_array)-too many tries?\n");
+ fprintf(E->trace.fpt,"%d %d %d\n",max_tries,number_of_tries,RAND_MAX);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+ random1=(1.0*rand())/(1.0*RAND_MAX);
+ random2=(1.0*rand())/(1.0*RAND_MAX);
+ random3=(1.0*rand())/(1.0*RAND_MAX);
+
+ x=dmin+random1*(dmax-dmin);
+ y=dmin+random2*(dmax-dmin);
+ z=dmin+random3*(dmax-dmin);
+
+/* first check if within shell */
+
+ rad=sqrt(x*x+y*y+z*z);
+
+ if (rad>=E->sx[j][3][E->lmesh.noz]) goto next_try;
+ if (rad<E->sx[j][3][1]) goto next_try;
+
+
+/* check if in cap */
+
+ ival=icheck_cap(E,ithiscap,x,y,z,rad);
+
+ if (ival!=1) goto next_try;
+
+/* Made it, so record tracer information */
+
+ cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
+
+ keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+
+ E->trace.itrac[j][0]++;
+ kk=E->trace.itrac[j][0];
+
+ E->trace.rtrac[j][0][kk]=theta;
+ E->trace.rtrac[j][1][kk]=phi;
+ E->trace.rtrac[j][2][kk]=rad;
+ E->trace.rtrac[j][3][kk]=x;
+ E->trace.rtrac[j][4][kk]=y;
+ E->trace.rtrac[j][5][kk]=z;
+
+
+ if (E->trace.itracer_type==1)
+ {
+ if (rad<=E->trace.z_interface) E->trace.etrac[j][0][kk]=1.0;
+ if (rad>E->trace.z_interface) E->trace.etrac[j][0][kk]=0.0;
+ }
+
+next_try:
+ ;
+ } /* end while */
+
+
+}/* end j */
+
+fprintf(stderr,"DONE Making Tracer Array (%d)\n",E->parallel.me);
+fflush(stderr);
+
+return;
+}
+
+/******** READ TRACER ARRAY *********************************************/
+/* */
+/* This function reads tracers from input file. */
+/* All processors read the same input file, then sort out which ones */
+/* belong. */
+
+void read_tracer_file(E)
+ struct All_variables *E;
+{
+
+char input_s[1000];
+
+int number_of_tracers;
+int kk;
+int icheck;
+int iestimate;
+int icushion;
+int ithiscap;
+int j;
+
+int icheck_cap();
+int icheck_processor_shell();
+int isum_tracers();
+void initialize_tracer_arrays();
+void keep_in_sphere();
+void sphere_to_cart();
+void cart_to_sphere();
+void expand_tracer_arrays();
+
+double x,y,z;
+double theta,phi,rad;
+double rdum1,rdum2,rdum3;
+
+FILE *fptracer;
+
+ fptracer=fopen(E->trace.tracer_file,"r");
+ fprintf(E->trace.fpt,"Opening %s\n",E->trace.tracer_file);
+
+ fgets(input_s,200,fptracer);
+ sscanf(input_s,"%d",&number_of_tracers);
+ fprintf(E->trace.fpt,"%d Tracers in file \n",number_of_tracers);
+
+/* initially size tracer arrays to number of tracers divided by processors */
+
+ icushion=100;
+
+ iestimate=number_of_tracers/E->parallel.nproc + icushion;
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ ithiscap=E->sphere.capid[j];
+
+ E->trace.itracsize[j]=iestimate;
+
+ initialize_tracer_arrays(E,j,iestimate);
+
+ for (kk=1;kk<=number_of_tracers;kk++)
+ {
+ fgets(input_s,200,fptracer);
+ sscanf(input_s,"%lf %lf %lf",&rdum1,&rdum2,&rdum3);
+
+ if (E->trace.icartesian_or_spherical_input==0)
+ {
+ x=rdum1;
+ y=rdum2;
+ z=rdum3;
+ cart_to_sphere(E,x,y,z,&theta,&phi,&rad);
+ }
+ if (E->trace.icartesian_or_spherical_input==1)
+ {
+ theta=rdum1;
+ phi=rdum2;
+ rad=rdum3;
+ sphere_to_cart(E,theta,phi,rad,&x,&y,&z);
+ }
+
+
+/* make sure theta, phi is in range, and radius is within bounds */
+
+ keep_in_sphere(E,&x,&y,&z,&theta,&phi,&rad);
+
+/* check whether tracer is within processor domain */
+
+ icheck=1;
+ if (E->parallel.nprocz>1) icheck=icheck_processor_shell(E,j,rad);
+ if (icheck!=1) goto next_tracer;
+
+ icheck=icheck_cap(E,ithiscap,x,y,z,rad);
+ if (icheck==0) goto next_tracer;
+
+/* if still here, tracer is in processor domain */
+
+
+ E->trace.itrac[j][0]++;
+
+ if (E->trace.itrac[j][0]>=(E->trace.itracsize[j]-5)) expand_tracer_arrays(E,j);
+
+ E->trace.rtrac[j][0][E->trace.itrac[j][0]]=theta;
+ E->trace.rtrac[j][1][E->trace.itrac[j][0]]=phi;
+ E->trace.rtrac[j][2][E->trace.itrac[j][0]]=rad;
+ E->trace.rtrac[j][3][E->trace.itrac[j][0]]=x;
+ E->trace.rtrac[j][4][E->trace.itrac[j][0]]=y;
+ E->trace.rtrac[j][5][E->trace.itrac[j][0]]=z;
+
+/* Assign composition for active tracers */
+
+ if (E->trace.itracer_type==1)
+ {
+ if (rad<=E->trace.z_interface) E->trace.etrac[j][0][kk]=1.0;
+ if (rad>E->trace.z_interface) E->trace.etrac[j][0][kk]=0.0;
+ }
+
+
+next_tracer:
+ ;
+ } /* end kk, number of tracers */
+
+ fprintf(E->trace.fpt,"Number of tracers in this cap is: %d\n",
+ E->trace.itrac[j][0]);
+
+ } /* end j */
+
+ fclose(fptracer);
+
+ icheck=isum_tracers(E);
+
+ if (icheck!=number_of_tracers)
+ {
+ fprintf(E->trace.fpt,"ERROR(read_tracer_file) - tracers != number in file\n");
+ fprintf(E->trace.fpt,"Tracers in system: %d\n", icheck);
+ fprintf(E->trace.fpt,"Tracers in file: %d\n", number_of_tracers);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+return;
+}
+
+
+/********** CART TO SPHERE ***********************/
+void cart_to_sphere(E,x,y,z,theta,phi,rad)
+ struct All_variables *E;
+ double x,y,z;
+ double *theta,*phi,*rad;
+{
+
+double temp;
+double myatan();
+
+ temp=x*x+y*y;
+
+ *rad=sqrt(temp+z*z);
+ *theta=atan2(sqrt(temp),z);
+ *phi=myatan(y,x);
+
+
+return;
+}
+
+/********** SPHERE TO CART ***********************/
+void sphere_to_cart(E,theta,phi,rad,x,y,z)
+ struct All_variables *E;
+ double theta,phi,rad;
+ double *x,*y,*z;
+{
+
+double sint,cost,sinf,cosf;
+double temp;
+
+ sint=sin(theta);
+ cost=cos(theta);
+ sinf=sin(phi);
+ cosf=cos(phi);
+
+ temp=rad*sint;
+
+ *x=temp*cosf;
+ *y=temp*sinf;
+ *z=rad*cost;
+
+
+return;
+}
+
+/********* ICHECK THAT PROCESSOR SHELL ********/
+/* */
+/* Checks whether a given radius is within */
+/* a given processors radial domain. */
+/* Returns 0 if not, 1 if so. */
+/* The domain is defined as including the bottom */
+/* radius, but excluding the top radius unless */
+/* we the processor domain is the one that */
+/* is at the surface (then both boundaries are */
+/* included). */
+
+
+int icheck_that_processor_shell(E,j,nprocessor,rad)
+ struct All_variables *E;
+ int j;
+ int nprocessor;
+ double rad;
+{
+
+
+double top_r,bottom_r;
+
+ top_r=E->trace.TOP_RAD[nprocessor][j];
+ bottom_r=E->trace.BOTTOM_RAD[nprocessor][j];
+
+
+/* First check bottom */
+
+ if (rad<bottom_r) return 0;
+
+/* Check top */
+
+
+ if (E->trace.ITOP[nprocessor][j]==1)
+ {
+ if (rad<=top_r) return 1;
+ }
+ else
+ {
+ if (rad<top_r) return 1;
+ }
+
+/* If here, means point is above processor */
+
+return 0;
+}
+
+
+/********** ICHECK PROCESSOR SHELL *************/
+/* returns -99 if rad is below correct shell */
+/* returns 0 if rad is above correct shell */
+/* returns 1 if correct shell */
+/* */
+/* Shell, here, refers to processor shell */
+/* */
+/* shell is defined as bottom boundary up to */
+/* and not including the top boundary unless */
+/* the shell in question is the top shell */
+
+int icheck_processor_shell(E,j,rad)
+ struct All_variables *E;
+ double rad;
+ int j;
+
+{
+
+static int top=0;
+static int bottom=0;
+static int been_here=0;
+double eps=0.0000001;
+
+
+ if (E->parallel.nprocz==1) return 1;
+
+ if (been_here==0)
+ {
+ if ((E->sx[j][3][1]-eps)<E->sphere.ri) bottom=1;
+ if ((E->sx[j][3][E->lmesh.noz]+eps)>E->sphere.ro) top=1;
+ been_here++;
+ }
+
+ if (top==1)
+ {
+ if (rad>=E->sx[j][3][1])
+ {
+ return 1;
+ }
+ return -99;
+ }
+
+ if (bottom==1)
+ {
+ if (rad<E->sx[j][3][E->lmesh.noz])
+ {
+ return 1;
+ }
+ return 0;
+ }
+
+ if ((rad>=E->sx[j][3][1])&&(rad<E->sx[j][3][E->lmesh.noz]))
+ {
+ return 1;
+ }
+ if (rad>=E->sx[j][3][E->lmesh.noz])
+ {
+ return 0;
+ }
+ if (rad<E->sx[j][3][1])
+ {
+ return -99;
+ }
+
+
+ fprintf(E->trace.fpt,"Should not be here\n");
+ fprintf(E->trace.fpt,"Error(check_shell) radius: %f\n",rad);
+ fflush(E->trace.fpt);
+ exit(10);
+
+}
+
+/******* ICHECK ELEMENT *************************************/
+/* */
+/* This function serves to determine if a point lies within */
+/* a given element */
+
+int icheck_element(E,j,nel,x,y,z,rad)
+ struct All_variables *E;
+ int nel,j;
+ double x,y,z,rad;
+{
+
+int icheck;
+int icheck_element_column();
+int icheck_shell();
+
+ icheck=icheck_shell(E,nel,rad);
+ if (icheck==0)
+ {
+ return 0;
+ }
+
+ icheck=icheck_element_column(E,j,nel,x,y,z,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 */
+
+int icheck_shell(E,nel,rad)
+ struct All_variables *E;
+ int nel;
+ double rad;
+{
+
+int ival;
+int ibottom_node;
+int itop_node;
+
+double bottom_rad;
+double 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];
+
+ ival=0;
+ if ((rad>=bottom_rad)&&(rad<top_rad)) ival=1;
+
+return ival;
+}
+
+/******** ICHECK ELEMENT COLUMN ****************************/
+/* */
+/* This function serves to determine if a point lies within */
+/* a given element's column */
+
+int icheck_element_column(E,j,nel,x,y,z,rad)
+ struct All_variables *E;
+ int nel,j;
+ double x,y,z,rad;
+{
+
+double test_point[4];
+double rnode[5][10];
+
+int ival;
+int kk;
+int node;
+int icheck_bounds();
+
+
+/* AKMA REMOVE FOR SPEED */
+
+ if ((nel<1) || (nel>E->lmesh.nel) )
+ {
+ fprintf(E->trace.fpt,"ERROR(icheck element column)-Weird nel: %d\n",nel);
+ fprintf(E->trace.fpt,"x: %f y: %f z: %f rad: %f\n",x,y,z,rad);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+
+ E->trace.istat_elements_checked++;
+
+/* surface coords of element nodes */
+
+ for (kk=1;kk<=4;kk++)
+ {
+
+ node=E->ien[j][nel].node[kk+4];
+
+ rnode[kk][1]=E->x[j][1][node];
+ rnode[kk][2]=E->x[j][2][node];
+ rnode[kk][3]=E->x[j][3][node];
+
+ rnode[kk][4]=E->sx[j][1][node];
+ rnode[kk][5]=E->sx[j][2][node];
+
+ rnode[kk][6]=E->trace.DSinCos[j][2][node]; /* cos(theta) */
+ rnode[kk][7]=E->trace.DSinCos[j][0][node]; /* sin(theta) */
+ rnode[kk][8]=E->trace.DSinCos[j][3][node]; /* cos(phi) */
+ rnode[kk][9]=E->trace.DSinCos[j][1][node]; /* sin(phi) */
+
+ }
+
+/* test_point - project to outer radius */
+
+ test_point[1]=x/rad;
+ test_point[2]=y/rad;
+ test_point[3]=z/rad;
+
+ ival=icheck_bounds(E,test_point,rnode[1],rnode[2],rnode[3],rnode[4]);
+
+
+return ival;
+}
+
+
+/********* ICHECK CAP ***************************************/
+/* */
+/* This function serves to determine if a point lies within */
+/* a given cap */
+/* */
+int icheck_cap(E,icap,x,y,z,rad)
+ struct All_variables *E;
+ int icap;
+ double x,y,z,rad;
+{
+
+double test_point[4];
+double rnode[5][10];
+
+int ival;
+int kk;
+int icheck_bounds();
+
+/* surface coords of cap nodes */
+
+
+ for (kk=1;kk<=4;kk++)
+ {
+
+ rnode[kk][1]=E->trace.XCAP[icap][kk];
+ rnode[kk][2]=E->trace.YCAP[icap][kk];
+ rnode[kk][3]=E->trace.ZCAP[icap][kk];
+ rnode[kk][4]=E->trace.THETA_CAP[icap][kk];
+ rnode[kk][5]=E->trace.PHI_CAP[icap][kk];
+ rnode[kk][6]=E->trace.COS_THETA[icap][kk];
+ rnode[kk][7]=E->trace.SIN_THETA[icap][kk];
+ rnode[kk][8]=E->trace.COS_PHI[icap][kk];
+ rnode[kk][9]=E->trace.SIN_PHI[icap][kk];
+ }
+
+
+/* test_point - project to outer radius */
+
+ test_point[1]=x/rad;
+ test_point[2]=y/rad;
+ test_point[3]=z/rad;
+
+ ival=icheck_bounds(E,test_point,rnode[1],rnode[2],rnode[3],rnode[4]);
+
+
+return ival;
+}
+
+/***** 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 */
+
+int icheck_bounds(E,test_point,rnode1,rnode2,rnode3,rnode4)
+ struct All_variables *E;
+ double *test_point;
+ double *rnode1;
+ double *rnode2;
+ double *rnode3;
+ double *rnode4;
+{
+
+int number_of_tries=0;
+int icheck;
+
+double v12[4];
+double v23[4];
+double v34[4];
+double v41[4];
+double v1p[4];
+double v2p[4];
+double v3p[4];
+double v4p[4];
+double cross1[4];
+double cross2[4];
+double cross3[4];
+double cross4[4];
+double rad1,rad2,rad3,rad4;
+double theta, phi;
+double tiny, eps;
+double x,y,z;
+
+double findradial();
+void makevector();
+void crossit();
+double myatan();
+
+/* make vectors from node to node */
+
+ makevector(v12,rnode2[1],rnode2[2],rnode2[3],rnode1[1],rnode1[2],rnode1[3]);
+ makevector(v23,rnode3[1],rnode3[2],rnode3[3],rnode2[1],rnode2[2],rnode2[3]);
+ makevector(v34,rnode4[1],rnode4[2],rnode4[3],rnode3[1],rnode3[2],rnode3[3]);
+ makevector(v41,rnode1[1],rnode1[2],rnode1[3],rnode4[1],rnode4[2],rnode4[3]);
+
+try_again:
+
+ number_of_tries++;
+
+/* make vectors from test point to node */
+
+ makevector(v1p,test_point[1],test_point[2],test_point[3],rnode1[1],rnode1[2],rnode1[3]);
+ makevector(v2p,test_point[1],test_point[2],test_point[3],rnode2[1],rnode2[2],rnode2[3]);
+ makevector(v3p,test_point[1],test_point[2],test_point[3],rnode3[1],rnode3[2],rnode3[3]);
+ makevector(v4p,test_point[1],test_point[2],test_point[3],rnode4[1],rnode4[2],rnode4[3]);
+
+/* Calculate cross products */
+
+ crossit(cross2,v12,v2p);
+ crossit(cross3,v23,v3p);
+ crossit(cross4,v34,v4p);
+ crossit(cross1,v41,v1p);
+
+/* Calculate radial component of cross products */
+
+ rad1=findradial(E,cross1,rnode1[6],rnode1[7],rnode1[8],rnode1[9]);
+ rad2=findradial(E,cross2,rnode2[6],rnode2[7],rnode2[8],rnode2[9]);
+ rad3=findradial(E,cross3,rnode3[6],rnode3[7],rnode3[8],rnode3[9]);
+ rad4=findradial(E,cross4,rnode4[6],rnode4[7],rnode4[8],rnode4[9]);
+
+/* 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[1],test_point[2],test_point[3]);
+ fprintf(E->trace.fpt,"Nodal points: 1: %f %f %f\n",rnode1[1],rnode1[2],rnode1[3]);
+ fprintf(E->trace.fpt,"Nodal points: 2: %f %f %f\n",rnode2[1],rnode2[2],rnode2[3]);
+ fprintf(E->trace.fpt,"Nodal points: 3: %f %f %f\n",rnode3[1],rnode3[2],rnode3[3]);
+ fprintf(E->trace.fpt,"Nodal points: 4: %f %f %f\n",rnode4[1],rnode4[2],rnode4[3]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ if (fabs(rad1)<=tiny||fabs(rad2)<=tiny||fabs(rad3)<=tiny||fabs(rad4)<=tiny)
+ {
+ x=test_point[1];
+ y=test_point[2];
+ z=test_point[3];
+ theta=myatan(sqrt(x*x+y*y),z);
+ phi=myatan(y,x);
+
+ if (theta<=M_PI/2.0)
+ {
+ theta=theta+eps;
+ }
+ else
+ {
+ theta=theta-eps;
+ }
+ phi=phi+eps;
+ x=sin(theta)*cos(phi);
+ y=sin(theta)*sin(phi);
+ z=cos(theta);
+ test_point[1]=x;
+ test_point[2]=y;
+ test_point[3]=z;
+
+ number_of_tries++;
+ goto try_again;
+
+ }
+
+ icheck=0;
+ if (rad1>0.0&&rad2>0.0&&rad3>0.0&&rad4>0.0) icheck=1;
+
+/*
+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);
+*/
+
+return icheck;
+
+}
+
+/****************************************************************************/
+/* FINDRADIAL */
+/* */
+/* This function finds the radial component of a Cartesian vector */
+
+
+double findradial(E,vec,cost,sint,cosf,sinf)
+ struct All_variables *E;
+ double *vec;
+ double cost,sint,cosf,sinf;
+{
+ double radialparti,radialpartj,radialpartk;
+ double radial;
+
+ radialparti=vec[1]*sint*cosf;
+ radialpartj=vec[2]*sint*sinf;
+ radialpartk=vec[3]*cost;
+
+ radial=radialparti+radialpartj+radialpartk;
+
+
+return radial;
+}
+
+
+/******************MAKEVECTOR*********************************************************/
+
+void makevector(vec,xf,yf,zf,x0,y0,z0)
+ double *vec;
+ double xf,yf,zf,x0,y0,z0;
+{
+
+ vec[1]=xf-x0;
+ vec[2]=yf-y0;
+ vec[3]=zf-z0;
+
+
+return;
+}
+
+/********************CROSSIT********************************************************/
+
+void crossit(cross,A,B)
+ double *cross;
+ double *A;
+ double *B;
+{
+
+ cross[1]=A[2]*B[3]-A[3]*B[2];
+ cross[2]=A[3]*B[1]-A[1]*B[3];
+ cross[3]=A[1]*B[2]-A[2]*B[1];
+
+
+return;
+}
+
+/*******************************************************************/
+/* Function SETUP DSINCOS */
+/* */
+/* open memory for and fill DSinCos */
+/* similar as E.SinCos but double instead of float */
+/* DSinCos[cap][0][node]=sint */
+/* DSinCos[cap][1][node]=sinf */
+/* DSinCos[cap][2][node]=cost */
+/* DSinCos[cap][3][node]=cosf */
+/* */
+
+void setup_dsincos(E)
+ struct All_variables *E;
+{
+
+ int j,kk;
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=0;kk<=3;kk++)
+ {
+ if ((E->trace.DSinCos[j][kk]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(setup_dsincos)-not enough memory(ib)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ for (kk=1;kk<=(E->lmesh.nno);kk++)
+ {
+ E->trace.DSinCos[j][0][kk]=sin(E->sx[j][1][kk]);
+ E->trace.DSinCos[j][1][kk]=sin(E->sx[j][2][kk]);
+ E->trace.DSinCos[j][2][kk]=cos(E->sx[j][1][kk]);
+ E->trace.DSinCos[j][3][kk]=cos(E->sx[j][2][kk]);
+ }
+
+ } /* end each cap */
+
+return;
+}
+
+/************ FIX RADIUS ********************************************/
+/* This function moves particles back in bounds if they left */
+/* during advection */
+
+void fix_radius(E,radius,theta,phi,x,y,z)
+ struct All_variables *E;
+ double *radius;
+ double *theta;
+ double *phi;
+ double *x;
+ double *y;
+ double *z;
+
+{
+
+double sint,cost,sinf,cosf,rad;
+
+static int been_here=0;
+static double max_radius;
+static double min_radius;
+
+ if (been_here==0)
+ {
+ max_radius=E->sphere.ro-E->trace.box_cushion;
+ min_radius=E->sphere.ri+E->trace.box_cushion;
+ been_here++;
+ }
+
+ if (*radius > max_radius)
+ {
+ *radius=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;
+ 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;
+ }
+
+
+return;
+}
+
+
+/******************************************************************/
+/* FIX PHI */
+/* */
+/* This function constrains the value of phi to be */
+/* between 0 and 2 PI */
+/* */
+
+void fix_phi(phi)
+ double *phi;
+
+{
+
+static double two_pi=2.0*M_PI;
+
+
+do_again:
+
+ if (*phi<0.0) *phi=*phi+two_pi;
+ else if (*phi>=two_pi) *phi=*phi-two_pi;
+
+ if ((*phi<0.0) || (*phi >=two_pi)) goto do_again;
+
+
+return;
+}
+
+/******************************************************************/
+/* FIX THETA */
+/* */
+/* This function constrains the value of theta to be */
+/* between 0 and PI */
+/* */
+void fix_theta(theta)
+ double *theta;
+{
+
+static double two_pi=2.0*M_PI;
+
+do_again2:
+
+ if (*theta<0.0) *theta=fabs(*theta);
+ else if (*theta>M_PI) *theta=two_pi-*theta;
+
+ if ((*theta<0.0) || (*theta >M_PI)) goto do_again2;
+
+
+
+return;
+}
+
+/********** IGET ELEMENT *****************************************/
+/* */
+/* This function returns the the real element for a given point. */
+/* Returns -99 in not in this cap. */
+/* iprevious_element, if known, is the last known element. If */
+/* it is not known, input a negative number. */
+
+int iget_element(E,j,iprevious_element,x,y,z,theta,phi,rad)
+ struct All_variables *E;
+ int j;
+ int iprevious_element;
+ double x,y,z;
+ double theta,phi,rad;
+{
+
+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;
+
+int iget_regel();
+int iquick_element_column_search();
+int icheck_cap();
+int icheck_regular_neighbors();
+int iget_radial_element();
+
+elx=E->lmesh.elx;
+ely=E->lmesh.ely;
+elz=E->lmesh.elz;
+
+
+ ntheta=0;
+ nphi=0;
+
+/* First check previous 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,theta,phi,&ntheta,&nphi);
+ if (iregel<=0)
+ {
+ return -99;
+ }
+
+
+/* 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,x,y,z,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,x,y,z,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,x,y,z,rad);
+ if (iel>0)
+ {
+ goto foundit;
+ }
+ }
+
+/* check if still in cap */
+
+ ival=icheck_cap(E,E->sphere.capid[j],x,y,z,rad);
+ if (ival==0)
+ {
+ return -99;
+ }
+
+/* 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],x,y,z,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,x,y,z,rad);
+ if (iel>0)
+ {
+ goto foundit;
+ }
+ }
+ }
+
+/* Decided to remove this part - not really needed and complicated */
+/*
+ else
+ {
+ iel=icheck_regular_neighbors(E,j,ntheta,nphi,x,y,z,theta,phi,rad);
+ if (iel>0)
+ {
+ goto foundit;
+ }
+ }
+*/
+ }
+
+/* As a last resort, check all element columns */
+
+ E->trace.istat1++;
+
+ iel=icheck_all_columns(E,j,x,y,z,rad);
+
+/*
+ fprintf(E->trace.fpt,"WARNING(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==1) exit(10);
+*/
+
+ if (E->trace.istat1%100==0)
+ {
+ fprintf(E->trace.fpt,"Checked all elements %d times already this turn\n",E->trace.istat1);
+ fprintf(stderr,"Checked all elements %d times already this turn\n",E->trace.istat1);
+ fflush(E->trace.fpt);
+ fflush(stderr);
+ }
+ if (iel>0)
+ {
+ goto foundit;
+ }
+
+
+/* if still here, there is a problem */
+
+ fprintf(E->trace.fpt,"Error(iget_element) - element not found\n");
+ fprintf(E->trace.fpt,"x,y,z,theta,phi,iregel %f %f %f %f %f %d\n",
+ x,y,z,theta,phi,iregel);
+ fflush(E->trace.fpt);
+ exit(10);
+
+foundit:
+
+ if (E->parallel.nprocz>1)
+ {
+ ival=icheck_processor_shell(E,j,rad);
+ if (ival!=1) return -99;
+ }
+
+/* find radial element */
+
+ ifinal_iel=iget_radial_element(E,j,iel,rad);
+
+return ifinal_iel;
+}
+
+
+/***** IGET RADIAL ELEMENT ***********************************/
+/* */
+/* This function returns the proper radial element, given */
+/* an element (iel) from the column. */
+
+int iget_radial_element(E,j,iel,rad)
+ struct All_variables *E;
+ int j,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;
+}
+
+/****** ICHECK REGULAR NEIGHBORS *****************************/
+/* */
+/* This function searches the regular element neighborhood. */
+
+/* This function is no longer used! */
+
+int icheck_regular_neighbors(E,j,ntheta,nphi,x,y,z,theta,phi,rad)
+ struct All_variables *E;
+ int j,ntheta,nphi;
+ double x,y,z;
+ double theta,phi,rad;
+{
+
+int new_ntheta,new_nphi;
+int kk,pp;
+int iregel;
+int ival;
+int imap[5];
+int ichoice;
+int irange;
+
+int iquick_element_column_search();
+
+fprintf(E->trace.fpt,"ERROR(icheck_regular_neighbors)-this subroutine is no longer used !\n");
+fflush(E->trace.fpt);
+exit(10);
+
+irange=2;
+
+ for (kk=-irange;kk<=irange;kk++)
+ {
+ for (pp=-irange;pp<=irange;pp++)
+ {
+ new_ntheta=ntheta+kk;
+ new_nphi=nphi+pp;
+ if ( (new_ntheta>0)&&(new_ntheta<=E->trace.numtheta[j])&&(new_nphi>0)&&(new_nphi<=E->trace.numphi[j]) )
+ {
+ iregel=new_ntheta+(new_nphi-1)*E->trace.numtheta[j];
+ if ((iregel>0) && (iregel<=E->trace.numregel[j]))
+ {
+ ival=iquick_element_column_search(E,j,iregel,new_ntheta,new_nphi,x,y,z,theta,phi,rad,imap,&ichoice);
+ if (ival>0) return ival;
+ }
+ }
+ }
+ }
+
+
+return -99;
+}
+
+
+
+
+
+/****** IQUICK ELEMENT SEARCH *****************************/
+/* */
+/* This function does a quick regular to real element */
+/* map check. Element number, if found, is returned. */
+/* Otherwise, -99 is returned. */
+/* Pointers to imap and ichoice are used because they may */
+/* prove to be convenient. */
+/* This routine is no longer used */
+
+int iquick_element_column_search(E,j,iregel,ntheta,nphi,x,y,z,theta,phi,rad,imap,ich)
+ struct All_variables *E;
+ int j,iregel;
+ int ntheta,nphi;
+ double x,y,z,theta,phi,rad;
+ int *imap;
+ int *ich;
+{
+
+int iregnode[5];
+int kk,pp;
+int nel,ival;
+int ichoice;
+int icount;
+int itemp1;
+int itemp2;
+
+int icheck_element_column();
+
+ fprintf(E->trace.fpt,"ERROR(iquick element)-this routine is no longer used!\n");
+ fflush(E->trace.fpt);
+ exit(10);
+
+/* REMOVE*/
+/*
+ ichoice=*ich;
+
+fprintf(E->trace.fpt,"AA: ichoice: %d\n",ichoice);
+fflush(E->trace.fpt);
+*/
+
+/* find regular nodes on regular element */
+
+/*
+ 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];
+*/
+
+ itemp1=iregel+nphi;
+ itemp2=itemp1+E->trace.numtheta[j];
+
+ iregnode[1]=itemp1-1;
+ iregnode[2]=itemp1;
+ iregnode[3]=itemp2+1;
+ iregnode[4]=itemp2;
+
+ for (kk=1;kk<=4;kk++)
+ {
+ if ((iregnode[kk]<1) || (iregnode[kk]>E->trace.numregnodes[j]) )
+ {
+ fprintf(E->trace.fpt,"ERROR(iquick)-weird regnode %d\n",iregnode[kk]);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+/* 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++;
+ imap[ichoice]=E->trace.regnodetoel[j][iregnode[kk]];
+
+
+next_corner:
+ ;
+ } /* end kk */
+
+ *ich=ichoice;
+
+/* statistical counter */
+
+ E->trace.istat_ichoice[j][ichoice]++;
+
+ if (ichoice==0) return -99;
+
+/* Here, no check is performed if all 4 corners */
+/* lie within a given element. */
+/* It may be possible (not sure) but unlikely */
+/* that the tracer is still not in that element */
+
+/* Decided to comment this out. */
+/* May not be valid for large regular grids. */
+/*
+*/
+/* AKMA */
+
+ if ((ichoice==1)&&(icount==4)) return imap[1];
+
+/* check others */
+
+ for (kk=1;kk<=ichoice;kk++)
+ {
+ nel=imap[kk];
+ ival=icheck_element_column(E,j,nel,x,y,z,rad);
+ if (ival>0) return nel;
+ }
+
+/* if still here, no element was found */
+
+return -99;
+}
+
+
+/*********** IGET REGEL ******************************************/
+/* */
+/* This function returns the regular element in which a point */
+/* exists. If not found, returns -99. */
+/* npi and ntheta are modified for later use */
+
+int iget_regel(E,j,theta,phi,ntheta,nphi)
+ struct All_variables *E;
+ int j;
+ double theta, 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 -99;
+ if (theta>E->trace.thetamax[j]) return -99;
+
+/* 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 -99;
+ if (iregel<1) return -99;
+
+return iregel;
+}
+
+
+/****** EXPAND TRACER ARRAYS *****************************************/
+
+void expand_tracer_arrays(E,j)
+ struct All_variables *E;
+ int j;
+
+{
+
+int inewsize;
+int kk;
+int icushion;
+
+/* expand rtrac and itrac by 20% */
+
+ icushion=100;
+
+ inewsize=E->trace.itracsize[j]+E->trace.itracsize[j]/5+icushion;
+
+ if ((E->trace.itrac[j]=(int *)realloc(E->trace.itrac[j],inewsize*sizeof(int)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (itracer)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ for (kk=0;kk<=((E->trace.number_of_advection_quantities)-1);kk++)
+ {
+ if ((E->trace.rtrac[j][kk]=(double *)realloc(E->trace.rtrac[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ for (kk=0;kk<=((E->trace.number_of_extra_quantities)-1);kk++)
+ {
+ if ((E->trace.etrac[j][kk]=(double *)realloc(E->trace.etrac[j][kk],inewsize*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(expand tracer arrays )-no memory 78 (%d)\n",kk);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+
+ fprintf(E->trace.fpt,"Expanding physical memory of itrac, rtrac, and etrac to %d from %d\n",
+ inewsize,E->trace.itracsize[j]);
+
+ E->trace.itracsize[j]=inewsize;
+
+return;
+}
+
+
+/**** DEBUG ***********************************************************/
+void debug(E,i)
+ struct All_variables *E;
+ int i;
+{
+
+ fprintf(E->trace.fpt,"HERE: %d\n",i);
+ fflush(E->trace.fpt);
+
+return;
+}
+
+/**** PDEBUG ***********************************************************/
+void pdebug(E,i)
+ struct All_variables *E;
+ int i;
+{
+
+ fprintf(E->trace.fpt,"HERE (Before Sync): %d\n",i);
+ fflush(E->trace.fpt);
+ parallel_process_sync();
+ fprintf(E->trace.fpt,"HERE (After Sync): %d\n",i);
+ fflush(E->trace.fpt);
+
+return;
+}
+
+/**** CDEBUG ***********************************************************/
+void cdebug(E,i)
+ struct All_variables *E;
+ int i;
+{
+
+ fprintf(stderr,"HERE (Before Sync): %d\n",i);
+ fflush(stderr);
+ parallel_process_sync();
+ fprintf(stderr,"HERE (After Sync): %d\n",i);
+ fflush(stderr);
+
+return;
+}
+
+/****************************************************************/
+/* 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 */
+/* UV[j][1][node]=u */
+/* UV[j][2][node]=v */
+
+void define_uv_space(E)
+ struct All_variables *E;
+
+{
+
+int kk,j;
+int midnode;
+int numnodes,node;
+
+double u,v,cosc,theta,phi;
+double theta_f,phi_f;
+
+ if (E->parallel.me==0) fprintf(stderr,"Setting up UV space\n");
+ fflush(stderr);
+
+ numnodes=E->lmesh.nno;
+
+/* open memory for uv space coords */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ for (kk=1;kk<=2;kk++)
+ {
+ if ((E->trace.UV[j][kk]=(double *)malloc((numnodes+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(define uv)-not enough memory(a)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+/* uv space requires a reference point */
+/* UV[j][1][0]=fixed theta */
+/* UV[j][2][0]=fixed phi */
+
+
+ midnode=numnodes/2;
+
+ E->trace.UV[j][1][0]=E->sx[j][1][midnode];
+ E->trace.UV[j][2][0]=E->sx[j][2][midnode];
+
+ theta_f=E->sx[j][1][midnode];
+ phi_f=E->sx[j][2][midnode];
+
+ E->trace.cos_theta_f=cos(theta_f);
+ E->trace.sin_theta_f=sin(theta_f);
+
+/* convert each nodal point to u and v */
+
+ for (node=1;node<=numnodes;node++)
+ {
+ theta=E->sx[j][1][node];
+ phi=E->sx[j][2][node];
+
+ cosc=cos(theta_f)*cos(theta)+sin(theta_f)*sin(theta)*
+ cos(phi-phi_f);
+ u=sin(theta)*sin(phi-phi_f)/cosc;
+ v=(sin(theta_f)*cos(theta)-cos(theta_f)*sin(theta)*cos(phi-phi_f))/cosc;
+
+ E->trace.UV[j][1][node]=u;
+ E->trace.UV[j][2][node]=v;
+
+ }
+
+
+ }/*end cap */
+
+return;
+}
+
+/***************************************************************/
+/* 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(E)
+ struct All_variables *E;
+{
+
+ int j,nelem,iwedge,kk;
+ int node;
+
+ double u[5],v[5];
+ double x1=0.0;
+ double x2=0.0;
+ double x3=0.0;
+ double y1=0.0;
+ double y2=0.0;
+ double y3=0.0;
+ double delta,a0,a1,a2;
+
+/* really only have to do this for each surface element, but */
+/* for simplicity, it is done for every element */
+
+if (E->parallel.me==0) fprintf(stderr," Determining Shape Coefficients\n");
+fflush(stderr);
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+/* first, allocate memory */
+
+ for(iwedge=1;iwedge<=2;iwedge++)
+ {
+ for (kk=1;kk<=9;kk++)
+ {
+ if ((E->trace.shape_coefs[j][iwedge][kk]=
+ (double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+ for (nelem=1;nelem<=E->lmesh.nel;nelem++)
+ {
+
+/* find u,v of local nodes at one radius */
+
+ for(kk=1;kk<=4;kk++)
+ {
+ node=E->IEN[E->mesh.levmax][j][nelem].node[kk];
+ u[kk]=E->trace.UV[j][1][node];
+ v[kk]=E->trace.UV[j][2][node];
+ }
+
+ for(iwedge=1;iwedge<=2;iwedge++)
+ {
+
+
+ if (iwedge==1)
+ {
+ x1=u[1];
+ x2=u[2];
+ x3=u[3];
+ y1=v[1];
+ y2=v[2];
+ y3=v[3];
+ }
+ if (iwedge==2)
+ {
+ x1=u[1];
+ x2=u[3];
+ x3=u[4];
+ y1=v[1];
+ y2=v[3];
+ y3=v[4];
+ }
+
+/* shape function 1 */
+
+ delta=(x3-x2)*(y1-y2)-(y2-y3)*(x2-x1);
+ a0=(x2*y3-x3*y2)/delta;
+ a1=(y2-y3)/delta;
+ a2=(x3-x2)/delta;
+
+ E->trace.shape_coefs[j][iwedge][1][nelem]=a0;
+ E->trace.shape_coefs[j][iwedge][2][nelem]=a1;
+ E->trace.shape_coefs[j][iwedge][3][nelem]=a2;
+
+/* shape function 2 */
+
+ delta=(x3-x1)*(y2-y1)-(y1-y3)*(x1-x2);
+ a0=(x1*y3-x3*y1)/delta;
+ a1=(y1-y3)/delta;
+ a2=(x3-x1)/delta;
+
+ E->trace.shape_coefs[j][iwedge][4][nelem]=a0;
+ E->trace.shape_coefs[j][iwedge][5][nelem]=a1;
+ E->trace.shape_coefs[j][iwedge][6][nelem]=a2;
+
+/* shape function 3 */
+
+ delta=(x1-x2)*(y3-y2)-(y2-y1)*(x2-x3);
+ a0=(x2*y1-x1*y2)/delta;
+ a1=(y2-y1)/delta;
+ a2=(x1-x2)/delta;
+
+ E->trace.shape_coefs[j][iwedge][7][nelem]=a0;
+ E->trace.shape_coefs[j][iwedge][8][nelem]=a1;
+ E->trace.shape_coefs[j][iwedge][9][nelem]=a2;
+
+
+ } /* end wedge */
+ } /* end elem */
+ } /* end cap */
+
+
+return;
+}
+
+/****************** GET CARTESIAN VELOCITY FIELD **************************************/
+/* */
+/* This function computes the cartesian velocity field from the spherical */
+/* velocity field computed from main Citcom code. */
+
+void get_cartesian_velocity_field(E)
+ struct All_variables *E;
+{
+
+int j,m,i;
+int kk;
+
+double sint,sinf,cost,cosf;
+double v_theta,v_phi,v_rad;
+double vx,vy,vz;
+
+static int been_here=0;
+
+if (been_here==0)
+{
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (kk=1;kk<=3;kk++)
+ {
+ if ((E->trace.V0_cart[j][kk]=(double *)malloc((E->lmesh.nno+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"ERROR(get_cartesian_velocity)-no memory 82hd7\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
+
+ been_here++;
+}
+
+
+for (m=1;m<=E->sphere.caps_per_proc;m++)
+{
+
+ for (i=1;i<=E->lmesh.nno;i++)
+ {
+ sint=E->trace.DSinCos[m][0][i];
+ sinf=E->trace.DSinCos[m][1][i];
+ cost=E->trace.DSinCos[m][2][i];
+ cosf=E->trace.DSinCos[m][3][i];
+
+ v_theta=E->sphere.cap[m].V[1][i];
+ v_phi=E->sphere.cap[m].V[2][i];
+ v_rad=E->sphere.cap[m].V[3][i];
+
+
+
+ vx=v_theta*cost*cosf-v_phi*sinf+v_rad*sint*cosf;
+ vy=v_theta*cost*sinf+v_phi*cosf+v_rad*sint*sinf;
+ vz=-v_theta*sint+v_rad*cost;
+
+ E->trace.V0_cart[m][1][i]=vx;
+ E->trace.V0_cart[m][2][i]=vy;
+ E->trace.V0_cart[m][3][i]=vz;
+ }
+}
+
+return;
+}
+
+/*********** KEEP IN SPHERE *********************************************/
+/* */
+/* This function makes sure the particle is within the sphere, and */
+/* phi and theta are within the proper degree range. */
+
+void keep_in_sphere(E,x,y,z,theta,phi,rad)
+ struct All_variable *E;
+ double *x;
+ double *y;
+ double *z;
+ double *theta;
+ double *phi;
+ double *rad;
+{
+
+void fix_radius();
+void fix_theta();
+void fix_phi();
+
+ fix_radius(E,rad,theta,phi,x,y,z);
+ fix_theta(theta);
+ fix_phi(phi);
+
+return;
+}
+
+/*********** CHECK SUM **************************************************/
+/* */
+/* This functions checks to make sure number of tracers is preserved */
+
+void check_sum(E)
+ struct All_variables *E;
+{
+
+int number,iold_number;
+
+static int been_here=0;
+
+int isum_tracers();
+
+ number=isum_tracers(E);
+
+ if (been_here==0)
+ {
+ E->trace.ilast_tracer_count=number;
+ fprintf(E->trace.fpt,"Sum of Tracers: %d\n",number);
+ been_here++;
+ return;
+ }
+
+ iold_number=E->trace.ilast_tracer_count;
+
+ if (number!=iold_number)
+ {
+ fprintf(E->trace.fpt,"ERROR(check_sum)-break in conservation %d %d\n",
+ number,iold_number);
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+
+ E->trace.ilast_tracer_count=number;
+
+return;
+}
+
+/************* ISUM TRACERS **********************************************/
+/* */
+/* This function uses MPI to sum all tracers and returns number of them. */
+
+int isum_tracers(E)
+ struct All_variables *E;
+{
+
+
+int imycount;
+int iallcount;
+int num;
+int j;
+
+ iallcount=0;
+
+imycount=0;
+for (j=1;j<=E->sphere.caps_per_proc;j++)
+{
+ imycount=imycount+E->trace.itrac[j][0];
+}
+
+ MPI_Allreduce(&imycount,&iallcount,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
+
+ num=iallcount;
+
+return num;
+}
+
+
+
+/* &&&&&&&&&&&&&&&&&&&& 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(E)
+ struct All_variables *E;
+
+{
+
+int kk,pp;
+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 vel_s[4];
+double vel_c[4];
+double my_theta0,my_phi0,my_rad0;
+double my_thetaf,my_phif,my_radf;
+double theta0,phi0,rad0;
+double thetaf,phif,radf;
+double x0_s[4],xf_s[4];
+double x0_c[4],xf_c[4];
+double vec[4];
+double runge_path_length,runge_time;
+double x0,y0,z0;
+double xf,yf,zf;
+double difference;
+double difperpath;
+
+void analytical_test_function();
+void predict_tracers();
+void correct_tracers();
+void analytical_runge_kutte();
+void sphere_to_cart();
+
+
+ fprintf(E->trace.fpt,"Starting Analytical Test\n");
+ if (E->parallel.me==0) fprintf(stderr,"Starting Analytical Test\n");
+ fflush(E->trace.fpt);
+ fflush(stderr);
+
+/* 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,theta,phi,rad,vel_s,vel_c);
+
+ E->sphere.cap[j].V[1][kk]=vel_s[1];
+ E->sphere.cap[j].V[2][kk]=vel_s[2];
+ E->sphere.cap[j].V[3][kk]=vel_s[3];
+ }
+ }
+
+ 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.itrac[j][0]>10)
+ {
+ fprintf(E->trace.fpt,"Warning(analytical)-too many tracers to print!\n");
+ fflush(E->trace.fpt);
+ if (E->trace.itracer_warnings==1) exit(10);
+ }
+ }
+
+/* print initial positions */
+
+ E->monitor.solution_cycles=0;
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (pp=1;pp<=E->trace.itrac[j][0];pp++)
+ {
+ theta=E->trace.rtrac[j][0][pp];
+ phi=E->trace.rtrac[j][1][pp];
+ rad=E->trace.rtrac[j][2][pp];
+
+ fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+
+ if (pp==1) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+
+ if (pp==1)
+ {
+ my_theta0=theta;
+ my_phi0=phi;
+ my_rad0=rad;
+ }
+ }
+ }
+
+/* advect tracers */
+
+ for (kk=1;kk<=nsteps;kk++)
+ {
+ E->monitor.solution_cycles=kk;
+
+ time=time+dt;
+
+ predict_tracers(E);
+ correct_tracers(E);
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (pp=1;pp<=E->trace.itrac[j][0];pp++)
+ {
+ theta=E->trace.rtrac[j][0][pp];
+ phi=E->trace.rtrac[j][1][pp];
+ rad=E->trace.rtrac[j][2][pp];
+
+ fprintf(E->trace.fpt,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+
+ if (pp==1) fprintf(stderr,"(%d) time: %f theta: %f phi: %f rad: %f\n",E->monitor.solution_cycles,time,theta,phi,rad);
+
+ if ((kk==nsteps) && (pp==1))
+ {
+ 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);
+ fflush(stderr);
+ parallel_process_sync();
+
+ 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.itrac[j][0];
+ }
+
+ MPI_Allreduce(&my_number,&number,1,MPI_INT,MPI_SUM,MPI_COMM_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);
+ fflush(stderr);
+ parallel_process_termination();
+ }
+
+
+/* communicate starting and final positions */
+
+ MPI_Allreduce(&my_theta0,&theta0,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+ MPI_Allreduce(&my_phi0,&phi0,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+ MPI_Allreduce(&my_rad0,&rad0,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+ MPI_Allreduce(&my_thetaf,&thetaf,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+ MPI_Allreduce(&my_phif,&phif,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+ MPI_Allreduce(&my_radf,&radf,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+
+ x0_s[1]=theta0;
+ x0_s[2]=phi0;
+ x0_s[3]=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 */
+
+ x0=x0_c[1];
+ y0=x0_c[2];
+ z0=x0_c[3];
+
+/* convert final citcom coords into cartesian */
+
+ sphere_to_cart(E,thetaf,phif,radf,&xf,&yf,&zf);
+
+ difference=sqrt((xf-xf_c[1])*(xf-xf_c[1])+(yf-xf_c[2])*(yf-xf_c[2])+(zf-xf_c[3])*(zf-xf_c[3]));
+
+ 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[1],xf_s[2],xf_s[3]);
+ 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[1],xf_s[2],xf_s[3]);
+ 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);
+ fflush(stderr);
+
+return;
+}
+
+/*************** ANALYTICAL RUNGE KUTTE ******************/
+/* */
+void analytical_runge_kutte(E,nsteps,dt,x0_s,x0_c,xf_s,xf_c,vec)
+ struct All_variables *E;
+ int nsteps;
+ double dt;
+ double *x0_c;
+ double *x0_s;
+ double *xf_c;
+ double *xf_s;
+ double *vec;
+
+{
+
+int kk;
+
+double x_0,y_0,z_0;
+double x_p,y_p,z_p;
+double x_c=0.0;
+double y_c=0.0;
+double z_c=0.0;
+double theta_0,phi_0,rad_0;
+double theta_p,phi_p,rad_p;
+double theta_c,phi_c,rad_c;
+double vel0_s[4],vel0_c[4];
+double velp_s[4],velp_c[4];
+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];
+
+ sphere_to_cart(E,theta_0,phi_0,rad_0,&x_0,&y_0,&z_0);
+
+/* fill initial cartesian vector to send back */
+
+ x0_c[1]=x_0;
+ x0_c[2]=y_0;
+ x0_c[3]=z_0;
+
+ time=0.0;
+ path=0.0;
+
+ for (kk=1;kk<=nsteps;kk++)
+ {
+
+/* get velocity at initial position */
+
+ analytical_test_function(E,theta_0,phi_0,rad_0,vel0_s,vel0_c);
+
+/* Find predicted midpoint position */
+
+ x_p=x_0+vel0_c[1]*dt*0.5;
+ y_p=y_0+vel0_c[2]*dt*0.5;
+ z_p=z_0+vel0_c[3]*dt*0.5;
+
+/* convert to spherical */
+
+ cart_to_sphere(E,x_p,y_p,z_p,&theta_p,&phi_p,&rad_p);
+
+/* get velocity at predicted midpoint position */
+
+ analytical_test_function(E,theta_p,phi_p,rad_p,velp_s,velp_c);
+
+/* Find corrected position using midpoint velocity */
+
+ x_c=x_0+velp_c[1]*dt;
+ y_c=y_0+velp_c[2]*dt;
+ z_c=z_0+velp_c[3]*dt;
+
+/* convert to spherical */
+
+ cart_to_sphere(E,x_c,y_c,z_c,&theta_c,&phi_c,&rad_c);
+
+/* compute path lenght */
+
+ dtpath=sqrt((x_c-x_0)*(x_c-x_0)+(y_c-y_0)*(y_c-y_0)+(z_c-z_0)*(z_c-z_0));
+ path=path+dtpath;
+
+ time=time+dt;
+
+ x_0=x_c;
+ y_0=y_c;
+ z_0=z_c;
+
+/* next time step */
+
+ }
+
+/* fill final spherical and cartesian vectors to send back */
+
+ xf_s[1]=theta_c;
+ xf_s[2]=phi_c;
+ xf_s[3]=rad_c;
+
+ xf_c[1]=x_c;
+ xf_c[2]=y_c;
+ xf_c[3]=z_c;
+
+ vec[1]=time;
+ vec[2]=path;
+
+return;
+}
+
+/**************** ANALYTICAL TEST FUNCTION ******************/
+/* */
+/* vel_s[1] => velocity in theta direction */
+/* vel_s[2] => velocity in phi direction */
+/* vel_s[3] => velocity in radial direction */
+/* */
+/* vel_c[1] => velocity in x direction */
+/* vel_c[2] => velocity in y direction */
+/* vel_c[3] => velocity in z direction */
+
+void analytical_test_function(E,theta,phi,rad,vel_s,vel_c)
+ struct All_variables *E;
+ double theta,phi,rad;
+ double *vel_s;
+ double *vel_c;
+
+{
+
+double sint,sinf,cost,cosf;
+double v_theta,v_phi,v_rad;
+double vx,vy,vz;
+
+/* This is where the function is given in spherical */
+
+ v_theta=50.0*rad*cos(phi);
+ v_phi=100.0*rad*sin(theta);
+ v_rad=25.0*rad;
+
+ vel_s[1]=v_theta;
+ vel_s[2]=v_phi;
+ vel_s[3]=v_rad;
+
+/* Convert the function into cartesian */
+
+ sint=sin(theta);
+ sinf=sin(phi);
+ cost=cos(theta);
+ cosf=cos(phi);
+
+ vx=v_theta*cost*cosf-v_phi*sinf+v_rad*sint*cosf;
+ vy=v_theta*cost*sinf+v_phi*cosf+v_rad*sint*sinf;
+ vz=-v_theta*sint+v_rad*cost;
+
+ vel_c[1]=vx;
+ vel_c[2]=vy;
+ vel_c[3]=vz;
+
+return;
+}
+
+/*************** TRACE FILTER ****************************************/
+/* */
+/* This function allows ad hoc reduction of entrainment. */
+/* This function should be used with caution. */
+/* itrace_filter==0 off */
+/* itrace_filter==1 simple filter */
+/* itrace_filter==2 relayering(see below) */
+/* */
+/* These methods are only to be used if compositions are 1.0 and 0.0 */
+
+void trace_filter(E)
+ struct All_variables *E;
+{
+
+int j;
+int kk;
+int numtracers;
+int nelem;
+int maxlevel=E->mesh.levmax;
+int ibottom_node;
+
+double element_composition;
+double element_volume;
+double dense_component_volume;
+double top_radius;
+double bottom_radius;
+double tracer_radius;
+double start_time;
+
+static int been_here=0;
+static double uplimit;
+static double downlimit;
+static double onethird=1.0/3.0;
+static int ibottom_elements;
+
+double CPU_time0();
+
+ if (!E->trace.itrace_filter) return;
+
+ if ((E->monitor.solution_cycles%E->trace.itrace_filter_every)!=0) return;
+
+ if (been_here==0)
+ {
+ if (E->trace.itrace_filter==1)
+ {
+ uplimit=E->trace.trace_filter_uplimit;
+ downlimit=E->trace.trace_filter_downlimit;
+ }
+ if (E->trace.itrace_filter==2)
+ {
+ ibottom_elements=12*E->lmesh.ely*E->lmesh.elx;
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ if ( (E->trace.trace_filter_toprad[j]=(double *) malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
+ {
+ fprintf(E->trace.fpt,"Error(trace_filter)-not enough memory!\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+
+ }
+ been_here++;
+ }
+
+ start_time=CPU_time0();
+
+/* Filter type 1 */
+
+ if (E->trace.itrace_filter==1)
+ {
+
+ fprintf(E->trace.fpt,"filtering tracers - type 1\n");
+ if (E->parallel.me==0) fprintf(stderr,"filtering tracers - type 1\n");
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ numtracers=E->trace.itrac[j][0];
+
+ for (kk=1;kk<=numtracers;kk++)
+ {
+ nelem=E->trace.itrac[j][kk];
+ element_composition=E->trace.comp_el[j][nelem];
+
+ if (element_composition>uplimit)
+ {
+ E->trace.etrac[j][0][kk]=1.0;
+ }
+ if (element_composition<downlimit)
+ {
+ E->trace.etrac[j][0][kk]=0.0;
+ }
+ }
+ }
+ }
+
+/* filter 2 seperates tracers with an element */
+
+ if (E->trace.itrace_filter==2)
+ {
+ fprintf(E->trace.fpt,"filtering tracers - type 2\n");
+ if (E->parallel.me==0) fprintf(stderr,"filtering tracers - type 2\n");
+
+/* scan all elements, determine volume of dense tracers */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+ for (nelem=1;nelem<=E->lmesh.nel;nelem++)
+ {
+ element_volume=E->ECO[maxlevel][j][nelem].area;
+ element_composition=E->trace.comp_el[j][nelem];
+ dense_component_volume=element_volume*element_composition;
+ ibottom_node=E->IEN[maxlevel][j][nelem].node[1];
+ bottom_radius=E->sx[j][3][ibottom_node];
+
+ top_radius=pow( (pow(bottom_radius,3.0)+dense_component_volume*ibottom_elements),(onethird));
+ E->trace.trace_filter_toprad[j][nelem]=top_radius;
+
+ }
+ }
+
+/* scan all tracers, adjust composition */
+
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ {
+
+ numtracers=E->trace.itrac[j][0];
+
+ for (kk=1;kk<=numtracers;kk++)
+ {
+ nelem=E->trace.itrac[j][kk];
+ top_radius=E->trace.trace_filter_toprad[j][nelem];
+ tracer_radius=E->trace.rtrac[j][2][kk];
+ if (tracer_radius>top_radius)
+ {
+ E->trace.etrac[j][0][kk]=0.0;
+ }
+ else if (tracer_radius<=top_radius)
+ {
+ E->trace.etrac[j][0][kk]=1.0;
+ }
+ }
+ }
+
+ }
+
+ fprintf(E->trace.fpt,"Filtering time: %f\n",CPU_time0()-start_time);
+ fflush(E->trace.fpt);
+
+return;
+}
+
+
+/******************* READ COMP ***********************************************/
+/* */
+/* This function is similar to read_temp. It is used to read the composition */
+/* from file for post-proceesing. */
+
+void read_comp(E)
+struct All_variables *E;
+ {
+ int i,ii,m,mm,ll;
+ char output_file[255],input_s[1000];
+
+ double g;
+ FILE *fp;
+
+
+
+ ii = E->monitor.solution_cycles;
+ sprintf(output_file,"%s.comp.%d.%d",E->control.old_P_file,E->parallel.me,ii);
+
+ if ((fp=fopen(output_file,"r"))==NULL)
+ {
+ fprintf(stderr,"ERROR(read_temp) - %s not found\n",output_file);
+ fflush(stderr);
+ exit(10);
+ }
+
+ fgets(input_s,1000,fp);
+ sscanf(input_s,"%d %d %f",&ll,&mm,&E->monitor.elapsed_time);
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ {
+ E->trace.comp_node[m]=(double *)malloc((E->lmesh.nno+1)*sizeof(double));
+
+ fgets(input_s,1000,fp);
+ sscanf(input_s,"%d %d",&ll,&mm);
+ for(i=1;i<=E->lmesh.nno;i++)
+ {
+ if (fgets(input_s,1000,fp)==NULL)
+ {
+ fprintf(stderr,"ERROR(read_comp) -data for node %d not found\n",i);
+ fflush(stderr);
+ exit(10);
+ }
+ sscanf(input_s,"%lf",&g);
+ E->trace.comp_node[m][i] = g;
+
+ }
+ }
+
+ fclose (fp);
+ return;
+}
+
+
+
+
+
+
+
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-03-01 19:36:33 UTC (rev 6144)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2007-03-01 21:02:00 UTC (rev 6145)
@@ -277,7 +277,150 @@
return((float)integral);
}
+/************ RETURN BULK VALUE_D *****************************************/
+/* */
+/* Same as Citcom function but allowing double instead of float. */
+/* I think when integer average =1, volume average is returned. */
+/* when integer average =0, integral is returned. */
+
+double return_bulk_value_d(E,Z,average)
+ struct All_variables *E;
+ double **Z;
+ int average;
+
+{
+ void get_global_shape_fn();
+ void float_global_operation();
+
+ int n,i,j,el,m;
+ double volume,integral,volume1,integral1;
+
+ const int vpts = vpoints[E->mesh.nsd];
+ const int ends = enodes[E->mesh.nsd];
+
+ volume1=0.0;
+ integral1=0.0;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for (el=1;el<=E->lmesh.nel;el++) {
+
+ for(j=1;j<=vpts;j++)
+ for(i=1;i<=ends;i++) {
+ n = E->ien[m][el].node[i];
+ volume1 += E->N.vpt[GNVINDEX(i,j)] * E->GDA[E->mesh.levmax][m][el].vpt[j-1];
+ integral1 += Z[m][n] * E->N.vpt[GNVINDEX(i,j)] * E->GDA[E->mesh.levmax][m][el].vpt[j-1];
+ }
+
+ }
+
+
+ MPI_Allreduce(&volume1 ,&volume ,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+ MPI_Allreduce(&integral1,&integral,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+
+ if(average && volume != 0.0)
+ integral /= volume;
+
+
+ return((double)integral);
+}
+
+/******* RETURN ELEMENTWISE HORIZ AVE ********************************/
+/* */
+/* This function is similar to return_horiz_ave in the citcom code */
+/* however here, elemental horizontal averages are given rather than */
+/* nodal averages. Also note, here is average per element */
+
+void return_elementwise_horiz_ave(E,X,H)
+ struct All_variables *E;
+ double **X, *H;
+{
+
+ int m,i,j,k,d,noz,noy,el,elz,elx,ely,nproc;
+ int sizeofH;
+ int elz2;
+ int *processors;
+ double *Have,*temp;
+
+ MPI_Comm world, horizon_p;
+ MPI_Group world_g, horizon_g;
+
+ sizeofH = (2*E->lmesh.elz+2)*sizeof(double);
+
+ processors = (int *)malloc((E->parallel.nprocxy+2)*sizeof(int));
+ Have = (double *)malloc(sizeofH);
+ temp = (double *)malloc(sizeofH);
+
+ noz = E->lmesh.noz;
+ noy = E->lmesh.noy;
+ elz = E->lmesh.elz;
+ elx = E->lmesh.elx;
+ ely = E->lmesh.ely;
+ elz2 = 2*elz;
+
+ for (i=0;i<=(elz*2+1);i++)
+ {
+ temp[i]=0.0;
+ }
+
+ for (i=1;i<=elz;i++)
+ {
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ {
+ for (k=1;k<=ely;k++)
+ {
+ for (j=1;j<=elx;j++)
+ {
+ el = i + (j-1)*elz + (k-1)*elx*elz;
+ temp[i] += X[m][el]*E->ECO[E->mesh.levmax][m][el].area;
+ temp[i+elz] += E->ECO[E->mesh.levmax][m][el].area;
+ }
+ }
+ }
+ }
+
+
+
+/* determine which processors should get the message from me for
+ computing the layer averages */
+
+ nproc = 0;
+ for (j=0;j<E->parallel.nprocxy;j++) {
+ d = E->parallel.me_loc[3] + j*E->parallel.nprocz;
+ processors[nproc] = d;
+ nproc ++;
+ }
+
+ if (nproc>0) {
+ world = MPI_COMM_WORLD;
+ MPI_Comm_group(world, &world_g);
+ MPI_Group_incl(world_g, nproc, processors, &horizon_g);
+ MPI_Comm_create(world, horizon_g, &horizon_p);
+
+ MPI_Allreduce(temp,Have,elz2+1,MPI_DOUBLE,MPI_SUM,horizon_p);
+
+ MPI_Comm_free (&horizon_p);
+ MPI_Group_free(&horizon_g);
+ MPI_Group_free(&world_g);
+
+ }
+ else
+ for (i=1;i<=elz2;i++)
+ Have[i] = temp[i];
+
+ for (i=1;i<=elz;i++) {
+ if(Have[i+elz] != 0.0)
+ H[i] = Have[i]/Have[i+elz];
+ }
+
+
+ free ((void *) Have);
+ free ((void *) temp);
+ free ((void *) processors);
+
+ return;
+}
+
/* ================================================== */
float find_max_horizontal(E,Tmax)
struct All_variables *E;
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2007-03-01 19:36:33 UTC (rev 6144)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2007-03-01 21:02:00 UTC (rev 6145)
@@ -115,6 +115,7 @@
Full_read_input_from_files.c \
Full_solver.c \
Full_sphere_related.c \
+ Full_tracer_advection.c \
Full_version_dependent.c \
Regional_boundary_conditions.c \
Regional_geometry_cartesian.c \
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-01 19:36:33 UTC (rev 6144)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-03-01 21:02:00 UTC (rev 6145)
@@ -685,7 +685,10 @@
struct SLICE slice;
struct Parallel parallel;
struct SPHERE sphere;
+ /* for regional tracer*/
struct Tracer Tracer;
+ /* for global tracer*/
+ struct TRACE trace;
struct Bdry boundary;
struct SBC sbc;
struct Output output;
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-01 19:36:33 UTC (rev 6144)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2007-03-01 21:02:00 UTC (rev 6145)
@@ -46,3 +46,128 @@
float *R_LOC_ELEM_T;
} r;
+
+
+struct TRACE{
+
+FILE *fpt;
+
+FILE *fp_error_fraction;
+FILE *fp_composition;
+
+char tracer_file[200];
+
+int number_of_advection_quantities;
+int number_of_extra_quantities;
+int number_of_tracer_quantities;
+int itracer_warnings;
+
+int itracing;
+int itracer_type;
+int ianalytical_tracer_test;
+int icompositional_rheology;
+int ibuoy_type;
+int itracer_restart;
+int icartesian_or_spherical_input;
+int itperel;
+int itracer_advection_scheme;
+int itracer_interpolation_scheme;
+int iwrite_tracers_every;
+int ireset_initial_composition;
+
+int itracsize[13];
+int *itrac[13];
+
+double buoyancy_ratio;
+double z_interface;
+double compositional_rheology_prefactor;
+double box_cushion;
+double *rtrac[13][100];
+double *etrac[13][100];
+double *rlater[13][100];
+/*
+double **rtrac[13];
+double **etrac[13];
+double **rlater[13];
+*/
+double *oldel[13];
+
+/* horizontally averaged */
+double *Have_C;
+double *Havel_tracers;
+
+/* tracer fileter */
+int itrace_filter;
+int itrace_filter_every;
+double trace_filter_uplimit;
+double trace_filter_downlimit;
+double *trace_filter_toprad[13];
+
+
+
+int ilater[13];
+int ilatersize[13];
+
+/* regular mesh parameters */
+int numtheta[13];
+int numphi[13];
+unsigned int numregel[13];
+unsigned int numregnodes[13];
+double deltheta[13];
+double delphi[13];
+double thetamax[13];
+double thetamin[13];
+double phimax[13];
+double phimin[13];
+int *regnodetoel[13];
+int *regtoel[13][5];
+
+
+/* statistical parameters */
+int istat_ichoice[13][5];
+int istat_isend;
+int istat_iempty;
+int istat1;
+int istat_elements_checked;
+int ilast_tracer_count;
+double stat_time_start_trace;
+double stat_time_end_trace;
+double stat_trace_time;
+double stat_stokes_time;
+
+/* Mesh information */
+double XCAP[13][5];
+double YCAP[13][5];
+double ZCAP[13][5];
+double THETA_CAP[13][5];
+double PHI_CAP[13][5];
+double RAD_CAP[13][5];
+double COS_THETA[13][5];
+double SIN_THETA[13][5];
+double COS_PHI[13][5];
+double SIN_PHI[13][5];
+double TOP_RAD[101][13];
+double BOTTOM_RAD[101][13];
+int ITOP[101][13];
+
+/* compositional information */
+int *ieltrac[13];
+double *celtrac[13];
+double *comp_el[13];
+double *comp_node[13];
+
+/* gnomonic shape functions */
+double *UV[13][3];
+double cos_theta_f;
+double sin_theta_f;
+double *shape_coefs[13][3][10];
+
+double *DSinCos[13][4];
+
+double *V0_cart[13][4];
+
+double initial_bulk_composition;
+double bulk_composition;
+double error_fraction;
+
+};
More information about the cig-commits
mailing list