[cig-commits] r5296 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Nov 15 17:52:26 PST 2006
Author: tan2
Date: 2006-11-15 17:52:25 -0800 (Wed, 15 Nov 2006)
New Revision: 5296
Modified:
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Full_obsolete.c
mc/3D/CitcomS/trunk/lib/Full_parallel_related.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Obsolete.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Regional_obsolete.c
mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Implemented geoid calculation, see issue64.
User visible changes include:
* "output_optional" can include string "geoid", which will enable geoid output
* new input parameters: radius, density_above, density_below
* removed input parameters: layerd, wdensity, llmax, nlong, nlati
* moved input parameter "output_ll_max" from Sphere.py to Output.py
* sync'd the defalut values of input parameters in Const.py and Instructions.c
TODO: h5output_geoid() not finished
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/* Functions to assemble the element k matrices and the element f vector.
@@ -77,7 +77,7 @@
const int nel=E->lmesh.nel;
const int lev=E->mesh.levmax;
- thermal_buoyancy(E,E->temp);
+ thermal_buoyancy(E,E->buoyancy);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -123,7 +123,7 @@
const int nel=E->lmesh.nel;
const int lev=E->mesh.levmax;
- thermal_buoyancy(E,E->temp);
+ thermal_buoyancy(E,E->buoyancy);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -765,7 +765,7 @@
for(p=0;p<n;p++) elt_f[p] = 0.0;
for(p=1;p<=ends;p++)
- force[p] = E->temp[m][E->ien[m][el].node[p]];
+ force[p] = E->buoyancy[m][E->ien[m][el].node[p]];
for(j=1;j<=vpts;j++) { /*compute force at each int point */
force_at_gs[j] = 0.0;
Modified: mc/3D/CitcomS/trunk/lib/Full_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Full_obsolete.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/*
@@ -591,70 +591,7 @@
}
-
/*************************************************************************/
-/* from Global_operations.c */
-/*************************************************************************/
-
-void sum_across_depth_sph1(E,sphc,sphs)
-struct All_variables *E;
-float *sphc,*sphs;
-{
- int jumpp,total,j,d;
-
- static float *sphcs,*temp;
- static int been_here=0;
- static int *processors,nproc;
-
- static MPI_Comm world, horizon_p;
- static MPI_Group world_g, horizon_g;
-
-if (been_here==0) {
- processors = (int *)malloc((E->parallel.nprocz+2)*sizeof(int));
- temp = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
- sphcs = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
-
- nproc = 0;
- for (j=0;j<E->parallel.nprocz;j++) {
- d =E->parallel.me_sph*E->parallel.nprocz+E->parallel.nprocz-1-j;
- processors[nproc] = d;
- nproc ++;
- }
-
- if (nproc>0) {
- world = E->parallel.world;
- MPI_Comm_group(world, &world_g);
- MPI_Group_incl(world_g, nproc, processors, &horizon_g);
- MPI_Comm_create(world, horizon_g, &horizon_p);
- }
-
- been_here++;
- }
-
- total = E->sphere.hindice*2+3;
- jumpp = E->sphere.hindice;
- for (j=0;j<E->sphere.hindice;j++) {
- sphcs[j] = sphc[j];
- sphcs[j+jumpp] = sphs[j];
- }
-
-
- if (nproc>0) {
-
- MPI_Allreduce(sphcs,temp,total,MPI_FLOAT,MPI_SUM,horizon_p);
-
- for (j=0;j<E->sphere.hindice;j++) {
- sphc[j] = temp[j];
- sphs[j] = temp[j+jumpp];
- }
-
- }
-
-return;
-}
-
-
-/*************************************************************************/
/* from Output.h */
/*************************************************************************/
@@ -843,7 +780,7 @@
fp1=fopen(output_file,"w");
fprintf(fp1,"lmaxx=%.4e lminx=%.4e for %s\n",maxx,minx,filen);
fprintf(fp1," ll mm cos sin \n");
- for (ll=0;ll<=E->sphere.output_llmax;ll++)
+ for (ll=0;ll<=E->output.llmax;ll++)
for(mm=0;mm<=ll;mm++) {
i = E->sphere.hindex[ll][mm];
fprintf(fp1,"%3d %3d %.4e %.4e \n",ll,mm,sphc[i],sphs[i]);
Modified: mc/3D/CitcomS/trunk/lib/Full_parallel_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_parallel_related.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Full_parallel_related.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -274,6 +274,9 @@
E->lmesh.snel = E->lmesh.elx*E->lmesh.ely;
+
+
+#if 0
i = cases[E->sphere.caps_per_proc];
E->parallel.nproc_sph[1] = incases3[i].xy[0];
@@ -297,8 +300,11 @@
E->sphere.lexs = E->sphere.lelx * E->parallel.me_loc_sph[1];
E->sphere.leys = E->sphere.lely * E->parallel.me_loc_sph[2];
+#endif
+
+
for(i=E->mesh.levmax;i>=E->mesh.levmin;i--) {
if (E->control.NMULTIGRID||E->control.EMULTIGRID) {
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -110,15 +110,20 @@
/* MPI_Barrier(E->parallel.world); */
}
+
+#if 0
E->sphere.elx = E->sphere.nox-1;
E->sphere.ely = E->sphere.noy-1;
E->sphere.snel = E->sphere.ely*E->sphere.elx;
E->sphere.nsf = E->sphere.noy*E->sphere.nox;
+#endif
+
+
/* Myr */
- E->data.scalet = (E->data.layer_km*1e3*E->data.layer_km*1e3/E->data.therm_diff)/(1.e6*365.25*24*3600);
+ E->data.scalet = (E->data.radius_km*1e3*E->data.radius_km*1e3/E->data.therm_diff)/(1.e6*365.25*24*3600);
/* cm/yr */
- E->data.scalev = (E->data.layer_km*1e3/E->data.therm_diff)/(100*365.25*24*3600);
+ E->data.scalev = (E->data.radius_km*1e3/E->data.therm_diff)/(100*365.25*24*3600);
E->data.timedir = E->control.Atemp / fabs(E->control.Atemp);
if(E->control.print_convergence && E->parallel.me==0)
Modified: mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <math.h>
@@ -238,15 +238,14 @@
counts =0;
if(E->parallel.me==0) fprintf(stderr,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
if(E->parallel.me==0) fprintf(E->fp,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
- valid = (residual < acc)?1:0;
- while (!valid) {
+ do {
residual=multi_grid(E,d0,F,acc,high_lev);
valid = (residual < acc)?1:0;
counts ++;
if(E->parallel.me==0) fprintf(stderr,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
if(E->parallel.me==0) fprintf(E->fp,"resi = %.6e for iter %d acc %.6e\n",residual,counts,acc);
- }
+ } while (!valid);
cycles = counts;
}
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <mpi.h>
@@ -325,6 +325,7 @@
temp = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
sphcs = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
+ /* pack */
jumpp = E->sphere.hindice;
total = E->sphere.hindice*2+3;
for (j=0;j<E->sphere.hindice;j++) {
@@ -332,8 +333,10 @@
sphcs[j+jumpp] = sphs[j];
}
+ /* sum across processors in horizontal direction */
MPI_Allreduce(sphcs,temp,total,MPI_FLOAT,MPI_SUM,E->parallel.horizontal_comm);
+ /* unpack */
for (j=0;j<E->sphere.hindice;j++) {
sphc[j] = temp[j];
sphs[j] = temp[j+jumpp];
@@ -613,4 +616,44 @@
temp1 = sqrt(temp2/temp1);
return (temp1);
- }
+}
+
+
+void sum_across_depth_sph1(E,sphc,sphs)
+ struct All_variables *E;
+ float *sphc,*sphs;
+{
+ int jumpp,total,j;
+
+ float *sphcs,*temp;
+
+ if (E->parallel.nprocz > 1) {
+ total = E->sphere.hindice*2+3;
+ temp = (float *) malloc(total*sizeof(float));
+ sphcs = (float *) malloc(total*sizeof(float));
+
+ /* pack sphc[] and sphs[] into sphcs[] */
+ jumpp = E->sphere.hindice;
+ for (j=0;j<E->sphere.hindice;j++) {
+ sphcs[j] = sphc[j];
+ sphcs[j+jumpp] = sphs[j];
+ }
+
+ /* sum across processors in z direction */
+ MPI_Allreduce(sphcs, temp, total, MPI_FLOAT, MPI_SUM,
+ E->parallel.vertical_comm);
+
+ /* unpack */
+ for (j=0;j<E->sphere.hindice;j++) {
+ sphc[j] = temp[j];
+ sphs[j] = temp[j+jumpp];
+ }
+
+ free(temp);
+ free(sphcs);
+ }
+
+
+ return;
+}
+
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -189,6 +189,7 @@
void set_convection_defaults();
void set_cg_defaults();
void set_mg_defaults();
+ float tmp;
int m=E->parallel.me;
/* first the problem type (defines subsequent behaviour) */
@@ -308,10 +309,7 @@
input_int("reset_startage",&(E->control.reset_startage),"0",m);
input_int("zero_elapsed_time",&(E->control.zero_elapsed_time),"0",m);
- input_int("ll_max",&(E->sphere.llmax),"1",m);
- input_int("nlong",&(E->sphere.noy),"1",m);
- input_int("nlati",&(E->sphere.nox),"1",m);
- input_int("output_ll_max",&(E->sphere.output_llmax),"1",m);
+ input_int("output_ll_max",&(E->output.llmax),"1",m);
input_int("topvbc",&(E->mesh.topvbc),"0",m);
input_int("botvbc",&(E->mesh.botvbc),"0",m);
@@ -363,22 +361,25 @@
/* data section */
input_float("Q0",&(E->control.Q0),"0.0",m);
- input_float("layerd",&(E->data.layer_km),"6371.0",m);
input_float("gravacc",&(E->data.grav_acc),"9.81",m);
- input_float("thermexp",&(E->data.therm_exp),"3.28e-5",m);
+ input_float("thermexp",&(E->data.therm_exp),"3.0e-5",m);
input_float("cp",&(E->data.Cp),"1200.0",m);
- input_float("thermdiff",&(E->data.therm_diff),"8.0e-7",m);
+ input_float("thermdiff",&(E->data.therm_diff),"1.0e-6",m);
input_float("density",&(E->data.density),"3340.0",m);
- input_float("wdensity",&(E->data.density_above),"1030.0",m);
+ input_float("density_above",&(E->data.density_above),"1030.0",m);
+ input_float("density_below",&(E->data.density_below),"6600.0",m);
input_float("refvisc",&(E->data.ref_viscosity),"1.0e21",m);
input_float("surftemp",&(E->data.surf_temp),"273.0",m);
+ input_float("radius",&tmp,"6371e3.0",m);
+ E->data.radius_km = tmp / 1e3;
+
E->data.therm_cond = E->data.therm_diff * E->data.density * E->data.Cp;
E->data.ref_temperature = E->control.Atemp * E->data.therm_diff
* E->data.ref_viscosity
/ (E->data.density * E->data.grav_acc * E->data.therm_exp)
- / (E->data.layer_km * E->data.layer_km * E->data.layer_km * 1e9);
+ / (E->data.radius_km * E->data.radius_km * E->data.radius_km * 1e9);
output_common_input(E);
h5input_params(E);
@@ -427,6 +428,7 @@
E->T[j] = (double *) malloc((nno+1)*sizeof(double));
E->NP[j] = (float *) malloc((nno+1)*sizeof(float));
E->edot[j] = (float *) malloc((nno+1)*sizeof(float));
+ E->buoyancy[j] = (double *) malloc((nno+1)*sizeof(double));
E->gstress[j] = (float *) malloc((6*nno+1)*sizeof(float));
E->stress[j] = (float *) malloc((12*nsf+1)*sizeof(float));
@@ -502,8 +504,9 @@
}
}
- for(i=0;i<=E->sphere.llmax;i++)
- E->sphere.hindex[i] = (int *) malloc((E->sphere.llmax+3)*sizeof(int));
+ for(i=0;i<=E->output.llmax;i++)
+ E->sphere.hindex[i] = (int *) malloc((E->output.llmax+3)
+ *sizeof(int));
for(i=E->mesh.gridmin;i<=E->mesh.gridmax;i++)
@@ -728,7 +731,6 @@
E->control.VBXbotval=0.0;
E->control.VBYbotval=0.0;
- E->data.layer_km = 2890.0; /* Earth, whole mantle defaults */
E->data.radius_km = 6370.0; /* Earth, whole mantle defaults */
E->data.grav_acc = 9.81;
E->data.therm_diff = 1.0e-6;
@@ -981,6 +983,7 @@
E->output.pressure = 0;
E->output.surf = 0;
E->output.botm = 0;
+ E->output.geoid = 0;
E->output.horiz_avg = 0;
while(1) {
@@ -1007,6 +1010,15 @@
E->output.surf = 1;
else if(strcmp(prev, "botm")==0)
E->output.botm = 1;
+ else if(strcmp(prev, "geoid")==0)
+ if (E->parallel.nprocxy != 12) {
+ fprintf(stderr, "Warning: geoid calculation only works in full version. Disabled\n");
+ }
+ else {
+ /* geoid calculation requires surface and CMB topo */
+ /* make sure the topos are available! */
+ E->output.geoid = 1;
+ }
else if(strcmp(prev, "horiz_avg")==0)
E->output.horiz_avg = 1;
else
Modified: mc/3D/CitcomS/trunk/lib/Obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Obsolete.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Obsolete.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -318,10 +318,391 @@
/* ========================================================== */
-/* */
+/* From Sphere_harmonics.c */
/* =========================================================== */
+
+/* This function construct sph harm tables on a regular grid */
+/* for inverse interpolation */
+
+static void compute_sphereh_int_table(E)
+ struct All_variables *E;
+{
+ int i,j;
+ double t,f;
+ double dth,dfi,sqrt_multis();
+
+ E->sphere.con = (double *)malloc((E->sphere.hindice+3)*sizeof(double));
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ E->sphere.con[E->sphere.hindex[ll][mm]] =
+ sqrt( (2.0-((mm==0)?1.0:0.0))*(2*ll+1)/(4.0*M_PI) )
+ *sqrt_multis(ll+mm,ll-mm); /* which is sqrt((ll-mm)!/(ll+mm)!) */
+ }
+
+ E->sphere.tablenplm = (double **) malloc((E->sphere.nox+1)
+ *sizeof(double*));
+ for (i=1;i<=E->sphere.nox;i++)
+ E->sphere.tablenplm[i]= (double *)malloc((E->sphere.hindice+3)
+ *sizeof(double));
+
+ E->sphere.tablencosf = (double **) malloc((E->sphere.noy+1)
+ *sizeof(double*));
+ E->sphere.tablensinf = (double **) malloc((E->sphere.noy+1)
+ *sizeof(double*));
+ for (i=1;i<=E->sphere.noy;i++) {
+ E->sphere.tablencosf[i]= (double *)malloc((E->output.llmax+3)
+ *sizeof(double));
+ E->sphere.tablensinf[i]= (double *)malloc((E->output.llmax+3)
+ *sizeof(double));
+ }
+
+ E->sphere.sx[1] = (double *) malloc((E->sphere.nsf+1)*sizeof(double));
+ E->sphere.sx[2] = (double *) malloc((E->sphere.nsf+1)*sizeof(double));
+
+ dth = M_PI/E->sphere.elx;
+ dfi = 2.0*M_PI/E->sphere.ely;
+
+ for (j=1;j<=E->sphere.noy;j++)
+ for (i=1;i<=E->sphere.nox;i++) {
+ node = i+(j-1)*E->sphere.nox;
+ E->sphere.sx[1][node] = dth*(i-1);
+ E->sphere.sx[2][node] = dfi*(j-1);
+ }
+
+ for (j=1;j<=E->sphere.nox;j++) {
+ t=E->sphere.sx[1][j];
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ p = E->sphere.hindex[ll][mm];
+ E->sphere.tablenplm[j][p] = modified_plgndr_a(ll,mm,t) ;
+ }
+ }
+ for (j=1;j<=E->sphere.noy;j++) {
+ node = 1+(j-1)*E->sphere.nox;
+ f=E->sphere.sx[2][node];
+ for (mm=0;mm<=E->output.llmax;mm++) {
+ E->sphere.tablencosf[j][mm] = cos( (double)(mm)*f );
+ E->sphere.tablensinf[j][mm] = sin( (double)(mm)*f );
+ }
+ }
+}
+
+
+/* ================================================
+ for a given node, this routine gives which cap and element
+ the node is in.
+ ================================================*/
+void construct_interp_net(E)
+ struct All_variables *E;
+{
+
+ void parallel_process_termination();
+ void parallel_process_sync();
+ int ii,jj,es,i,j,m,el,node;
+ int locate_cap(),locate_element();
+ double x[4],t,f;
+
+ const int ends=4;
+
+ for (i=1;i<=E->sphere.nox;i++)
+ for (j=1;j<=E->sphere.noy;j++) {
+ node = i+(j-1)*E->sphere.nox;
+ E->sphere.int_cap[node]=0;
+ E->sphere.int_ele[node]=0;
+ }
+
+
+ for (i=1;i<=E->sphere.nox;i++)
+ for (j=1;j<=E->sphere.noy;j++) {
+ node = i+(j-1)*E->sphere.nox;
+
+ /* first find which cap this node (i,j) is in */
+ t = E->sphere.sx[1][node];
+ f = E->sphere.sx[2][node];
+
+ x[1] = sin(t)*cos(f); /* radius does not matter */
+ x[2] = sin(t)*sin(f);
+ x[3] = cos(t);
+
+
+ fprintf(E->fp,"mmm0=%d\n",node);
+
+ m = locate_cap(E,x);
+
+ fprintf(E->fp,"mmm=%d\n",m);
+
+ if (m>0) {
+ el = locate_element(E,m,x,node); /* bottom element */
+
+ if (el<=0) {
+ fprintf(stderr,"!!! Processor %d cannot find the right element in cap %d\n",E->parallel.me,m);
+ parallel_process_termination();
+ }
+
+ E->sphere.int_cap[node]=m;
+ E->sphere.int_ele[node]=el;
+
+ }
+ } /* end for i and j */
+
+ parallel_process_sync(E);
+
+ return;
+}
+
+/* ================================================
+ locate the cap for node (i,j)
+ ================================================*/
+
+int locate_cap(E,x)
+ struct All_variables *E;
+ double x[4];
+{
+
+ int ia[5],i,m,mm;
+ double xx[4],angle[5],angle1[5];
+ double get_angle();
+ double area1,rr;
+ const double e_7=1.e-7;
+ static int been_here=0;
+
+ mm = 0;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ ia[1] = 1;
+ ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
+ ia[3] = E->lmesh.nno-E->lmesh.noz+1;
+ ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
+
+ for (i=1;i<=4;i++) {
+ xx[1] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
+ xx[2] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
+ xx[3] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
+ angle[i] = get_angle(x,xx); /* get angle between (i,j) and other four*/
+ angle1[i]=E->sphere.angle[m][i];
+ }
+
+ area1 = area_of_sphere_triag(angle[1],angle[2],angle1[1])
+ + area_of_sphere_triag(angle[2],angle[3],angle1[2])
+ + area_of_sphere_triag(angle[3],angle[4],angle1[3])
+ + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
+
+ if ( fabs ((area1-E->sphere.area[m])/E->sphere.area[m]) <e_7 ) {
+ mm = m;
+ return (mm);
+ }
+ }
+
+ return (mm);
+}
+
+/* ================================================
+ locate the element containing the node (i,j) with coord x.
+ The radius is assumed to be 1 in computing the areas.
+ NOTE: The returned element el is for the bottom layer.
+ ================================================*/
+
+int locate_element(E,m,x,ne)
+ struct All_variables *E;
+ double x[4];
+int m,ne;
+{
+
+ int el_temp,el,es,el_located,level,lev,lev_plus,el_plus,es_plus,i,j,node;
+ double area,area1,areamin;
+ double area_of_5points();
+ const double e_7=1.e-7;
+ const double e_6=1.e6;
+
+ el_located = 0;
+
+
+ level=E->mesh.levmin;
+ for (es=1;es<=E->lmesh.SNEL[level];es++) {
+
+ el = (es-1)*E->lmesh.ELZ[level]+1;
+ area1 = area_of_5points (E,level,m,el,x,ne);
+ area = E->sphere.area1[level][m][es];
+
+ if(fabs ((area1-area)/area) <e_7 ) {
+ for (lev=E->mesh.levmin;lev<E->mesh.levmax;lev++) {
+ lev_plus = lev + 1;
+ j=1;
+ areamin = e_6;
+ do {
+ el_plus = E->EL[lev][m][el].sub[j];
+
+ es_plus = (el_plus-1)/E->lmesh.ELZ[lev_plus]+1;
+
+ area1 = area_of_5points(E,lev_plus,m,el_plus,x,ne);
+ area = E->sphere.area1[lev_plus][m][es_plus];
+
+ if(fabs(area1-area)<areamin) {
+ areamin=fabs(area1-area);
+ el_temp = el_plus;
+ }
+ j++;
+ } while (j<5 && fabs((area1-area)/area) > e_7);
+ el = el_plus;
+ /* if exit with ..>e_7, pick the best one*/
+ if (fabs((area1-area)/area) > e_7) el = el_temp;
+ } /* end for loop lev */
+ el_located = el;
+ } /* end for if */
+
+ if(el_located) break;
+ } /* end for es at the coarsest level */
+
+ return (el_located);
+}
+
+/* ===============================================================
+ interpolate nodal T's within cap m and element el onto node with
+ coordinate x[3] which is derived from a regular mesh and within
+ the element el. NOTE the radius of x[3] is the inner radius.
+ =============================================================== */
+
+float sphere_interpolate_point(E,T,m,el,x,ne)
+ struct All_variables *E;
+ float **T;
+ double x[4];
+int m,el,ne;
+{
+ double to,fo,y[4],yy[4][5],dxdy[4][4];
+ double a1,b1,c1,d1,a2,b2,c2,d2,a,b,c,xx1,yy1,y1,y2;
+ float ta,t[5];
+ int es,snode,i,j,node;
+
+ const int oned=4;
+ const double e_7=1.e-7;
+ const double four=4.0;
+ const double two=2.0;
+ const double one=1.0;
+ const double pt25=0.25;
+
+ /* first rotate the coord such that the center of element is
+ the pole */
+
+ es = (el-1)/E->lmesh.elz+1;
+
+ to = E->eco[m][el].centre[1];
+ fo = E->eco[m][el].centre[2];
+
+ dxdy[1][1] = cos(to)*cos(fo);
+ dxdy[1][2] = cos(to)*sin(fo);
+ dxdy[1][3] = -sin(to);
+ dxdy[2][1] = -sin(fo);
+ dxdy[2][2] = cos(fo);
+ dxdy[2][3] = 0.0;
+ dxdy[3][1] = sin(to)*cos(fo);
+ dxdy[3][2] = sin(to)*sin(fo);
+ dxdy[3][3] = cos(to);
+
+ for(i=1;i<=oned;i++) { /* nodes */
+ node = E->ien[m][el].node[i];
+ snode = E->sien[m][es].node[i];
+ t[i] = T[m][snode];
+ for (j=1;j<=E->mesh.nsd;j++)
+ yy[j][i] = E->x[m][1][node]*dxdy[j][1]
+ + E->x[m][2][node]*dxdy[j][2]
+ + E->x[m][3][node]*dxdy[j][3];
+ }
+
+ for (j=1;j<=E->mesh.nsd;j++)
+ y[j] = x[1]*dxdy[j][1] + x[2]*dxdy[j][2] + x[3]*dxdy[j][3];
+
+ /* then for node y, determine its coordinates xx1,yy1
+ in the parental element in the isoparametric element system*/
+
+ a1 = yy[1][1] + yy[1][2] + yy[1][3] + yy[1][4];
+ b1 = yy[1][3] + yy[1][2] - yy[1][1] - yy[1][4];
+ c1 = yy[1][3] + yy[1][1] - yy[1][2] - yy[1][4];
+ d1 = yy[1][3] + yy[1][4] - yy[1][1] - yy[1][2];
+ a2 = yy[2][1] + yy[2][2] + yy[2][3] + yy[2][4];
+ b2 = yy[2][3] + yy[2][2] - yy[2][1] - yy[2][4];
+ c2 = yy[2][3] + yy[2][1] - yy[2][2] - yy[2][4];
+ d2 = yy[2][3] + yy[2][4] - yy[2][1] - yy[2][2];
+
+ a = d2*c1;
+ b = a2*c1+b1*d2-d1*c2-d1*b2-four*c1*y[2];
+ c=four*c2*y[1]-c2*a1-a1*b2+four*b2*y[1]-four*b1*y[2]+a2*b1;
+
+ if (fabs(a)<e_7) {
+ yy1 = -c/b;
+ xx1 = (four*y[1]-a1-d1*yy1)/(b1+c1*yy1);
+ }
+ else {
+ y1= (-b+sqrt(b*b-four*a*c))/(two*a);
+ y2= (-b-sqrt(b*b-four*a*c))/(two*a);
+ if (fabs(y1)>fabs(y2))
+ yy1 = y2;
+ else
+ yy1 = y1;
+ xx1 = (four*y[1]-a1-d1*yy1)/(b1+c1*yy1);
+ }
+
+ /* now we can calculate T at x[4] using shape function */
+
+ ta = ((one-xx1)*(one-yy1)*t[1]+(one+xx1)*(one-yy1)*t[2]+
+ (one+xx1)*(one+yy1)*t[3]+(one-xx1)*(one+yy1)*t[4])*pt25;
+
+ /*if(fabs(xx1)>1.5 || fabs(yy1)>1.5)fprintf(E->fp_out,"ME= %d %d %d %g %g %g %g %g %g %g\n",ne,m,es,t[1],t[2],t[3],t[4],ta,xx1,yy1);
+ */
+ return (ta);
+}
+
+/* ===================================================================
+ do the interpolation on sphere for data T, which is needed for both
+ spherical harmonic expansion and graphics
+ =================================================================== */
+
+void sphere_interpolate(E,T,TG)
+ struct All_variables *E;
+ float **T,*TG;
+{
+
+ float sphere_interpolate_point();
+ void gather_TG_to_me0();
+ void parallel_process_termination();
+
+ int ii,jj,es,i,j,m,el,node;
+ double x[4],t,f;
+
+ const int ends=4;
+
+ TG[0] = 0.0;
+ for (i=1;i<=E->sphere.nox;i++)
+ for (j=1;j<=E->sphere.noy;j++) {
+ node = i+(j-1)*E->sphere.nox;
+ TG[node] = 0.0;
+ /* first find which cap this node (i,j) is in */
+
+ m = E->sphere.int_cap[node];
+ el = E->sphere.int_ele[node];
+
+ if (m>0 && el>0) {
+ t = E->sphere.sx[1][node];
+ f = E->sphere.sx[2][node];
+
+ x[1] = E->sx[1][3][1]*sin(t)*cos(f);
+ x[2] = E->sx[1][3][1]*sin(t)*sin(f);
+ x[3] = E->sx[1][3][1]*cos(t);
+
+ TG[node] = sphere_interpolate_point(E,T,m,el,x,node);
+
+ }
+
+ } /* end for i and j */
+
+ gather_TG_to_me0(E,TG);
+
+ return;
+}
+
+
+
/* version */
/* $Id$ */
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -42,6 +42,7 @@
void output_visc_prepare(struct All_variables *, float **);
void output_visc(struct All_variables *, int);
void output_surf_botm(struct All_variables *, int);
+void output_geoid(struct All_variables *, int);
void output_stress(struct All_variables *, int);
void output_horiz_avg(struct All_variables *, int);
void output_tracer(struct All_variables *, int);
@@ -50,7 +51,7 @@
extern void parallel_process_termination();
extern void heat_flux(struct All_variables *);
extern void get_STD_topo(struct All_variables *, float**, float**,
- float**, float**, int);
+ float**, float**, int);
/**********************************************************************/
@@ -83,6 +84,10 @@
output_tracer(E, cycles);
/* optiotnal output below */
+ /* compute and output geoid (in spherical harmonics coeff) */
+ if (E->output.geoid == 1)
+ output_geoid(E, cycles);
+
if (E->output.stress == 1)
output_stress(E, cycles);
@@ -144,7 +149,7 @@
int lev = E->mesh.levmax;
sprintf(output_file,"%s.visc.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp1 = output_open(output_file);
@@ -167,7 +172,7 @@
FILE *fp1;
sprintf(output_file,"%s.velo.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp1 = output_open(output_file);
fprintf(fp1,"%d %d %.5e\n",cycles,E->lmesh.nno,E->monitor.elapsed_time);
@@ -197,21 +202,21 @@
if (E->output.surf && (E->parallel.me_loc[3]==E->parallel.nprocz-1)) {
sprintf(output_file,"%s.surf.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp2 = output_open(output_file);
for(j=1;j<=E->sphere.caps_per_proc;j++) {
- /* choose either STD topo or pseudo-free-surf topo */
- if(E->control.pseudo_free_surf)
- topo = E->slice.freesurf[j];
- else
- topo = E->slice.tpg[j];
+ /* choose either STD topo or pseudo-free-surf topo */
+ if(E->control.pseudo_free_surf)
+ topo = E->slice.freesurf[j];
+ else
+ topo = E->slice.tpg[j];
- fprintf(fp2,"%3d %7d\n",j,E->lmesh.nsf);
- for(i=1;i<=E->lmesh.nsf;i++) {
- s = i*E->lmesh.noz;
- fprintf(fp2,"%.4e %.4e %.4e %.4e\n",topo[i],E->slice.shflux[j][i],E->sphere.cap[j].V[1][s],E->sphere.cap[j].V[2][s]);
- }
+ fprintf(fp2,"%3d %7d\n",j,E->lmesh.nsf);
+ for(i=1;i<=E->lmesh.nsf;i++) {
+ s = i*E->lmesh.noz;
+ fprintf(fp2,"%.4e %.4e %.4e %.4e\n",topo[i],E->slice.shflux[j][i],E->sphere.cap[j].V[1][s],E->sphere.cap[j].V[2][s]);
+ }
}
fclose(fp2);
}
@@ -219,13 +224,13 @@
if (E->output.botm && (E->parallel.me_loc[3]==0)) {
sprintf(output_file,"%s.botm.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp2 = output_open(output_file);
for(j=1;j<=E->sphere.caps_per_proc;j++) {
fprintf(fp2,"%3d %7d\n",j,E->lmesh.nsf);
for(i=1;i<=E->lmesh.nsf;i++) {
- s = (i-1)*E->lmesh.noz + 1;
+ s = (i-1)*E->lmesh.noz + 1;
fprintf(fp2,"%.4e %.4e %.4e %.4e\n",E->slice.tpgb[j][i],E->slice.bhflux[j][i],E->sphere.cap[j].V[1][s],E->sphere.cap[j].V[2][s]);
}
}
@@ -236,6 +241,41 @@
}
+void output_geoid(struct All_variables *E, int cycles)
+{
+ void compute_geoid();
+ int ll, mm, p;
+ char output_file[255];
+ FILE *fp1;
+
+ compute_geoid(E, E->sphere.harm_geoid,
+ E->sphere.harm_tpgt, E->sphere.harm_tpgb);
+
+ if (E->parallel.me == 0) {
+ sprintf(output_file, "%s.geoid.%d.%d", E->control.data_file,
+ E->parallel.me, cycles);
+ fp1 = output_open(output_file);
+
+ /* write headers */
+ fprintf(fp1, "%d %d %.5e\n", cycles, E->output.llmax,
+ E->monitor.elapsed_time);
+
+ /* write sph harm coeff of geoid and topos */
+ for (ll=0; ll<=E->output.llmax; ll++)
+ for(mm=0; mm<=ll; mm++) {
+ p = E->sphere.hindex[ll][mm];
+ fprintf(fp1,"%d %d %.4e %.4e %.4e %.4e\n",
+ ll, mm,
+ E->sphere.harm_geoid[0][p], E->sphere.harm_geoid[1][p],
+ E->sphere.harm_tpgt[0][p], E->sphere.harm_tpgt[1][p]);
+ }
+
+ fclose(fp1);
+ }
+}
+
+
+
void output_stress(struct All_variables *E, int cycles)
{
int m, node;
@@ -243,7 +283,7 @@
FILE *fp1;
sprintf(output_file,"%s.stress.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp1 = output_open(output_file);
fprintf(fp1,"%d %d %.5e\n",cycles,E->lmesh.nno,E->monitor.elapsed_time);
@@ -252,12 +292,12 @@
fprintf(fp1,"%3d %7d\n",m,E->lmesh.nno);
for (node=1;node<=E->lmesh.nno;node++)
fprintf(fp1, "%.4e %.4e %.4e %.4e %.4e %.4e\n",
- E->gstress[m][(node-1)*6+1],
- E->gstress[m][(node-1)*6+2],
- E->gstress[m][(node-1)*6+3],
- E->gstress[m][(node-1)*6+4],
- E->gstress[m][(node-1)*6+5],
- E->gstress[m][(node-1)*6+6]);
+ E->gstress[m][(node-1)*6+1],
+ E->gstress[m][(node-1)*6+2],
+ E->gstress[m][(node-1)*6+3],
+ E->gstress[m][(node-1)*6+4],
+ E->gstress[m][(node-1)*6+5],
+ E->gstress[m][(node-1)*6+6]);
}
fclose(fp1);
}
@@ -279,7 +319,7 @@
if (E->parallel.me<E->parallel.nprocz) {
sprintf(output_file,"%s.horiz_avg.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp1=fopen(output_file,"w");
for(j=1;j<=E->lmesh.noz;j++) {
fprintf(fp1,"%.4e %.4e %.4e %.4e\n",E->sx[1][3][j],E->Have.T[j],E->Have.V[1][j],E->Have.V[2][j]);
@@ -319,7 +359,7 @@
FILE *fp1;
sprintf(output_file,"%s.pressure.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp1 = output_open(output_file);
fprintf(fp1,"%d %d %.5e\n",cycles,E->lmesh.nno,E->monitor.elapsed_time);
@@ -344,7 +384,7 @@
FILE *fp1;
sprintf(output_file,"%s.tracer.%d.%d", E->control.data_file,
- E->parallel.me, cycles);
+ E->parallel.me, cycles);
fp1 = output_open(output_file);
fprintf(fp1,"%.5e\n",E->monitor.elapsed_time);
@@ -367,11 +407,11 @@
if (E->parallel.me == 0) {
fprintf(E->fptime,"%d %.4e %.4e %.4e %.4e\n",
- cycles,
- E->monitor.elapsed_time,
- E->advection.timestep,
- current_time - E->monitor.cpu_time_at_start,
- current_time - E->monitor.cpu_time_at_last_cycle);
+ cycles,
+ E->monitor.elapsed_time,
+ E->advection.timestep,
+ current_time - E->monitor.cpu_time_at_start,
+ current_time - E->monitor.cpu_time_at_last_cycle);
fflush(E->fptime);
}
@@ -380,3 +420,4 @@
return;
}
+
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -87,6 +87,7 @@
void h5output_stress(struct All_variables *, int);
void h5output_tracer(struct All_variables *, int);
void h5output_surf_botm(struct All_variables *, int);
+void h5output_geoid(struct All_variables *, int);
void h5output_horiz_avg(struct All_variables *, int);
void h5output_time(struct All_variables *, int);
@@ -131,6 +132,9 @@
h5output_tracer(E, cycles);
/* optional output below */
+ if(E->output.geoid == 1)
+ h5output_geoid(E, cycles);
+
if(E->output.stress == 1)
h5output_stress(E, cycles);
@@ -1583,6 +1587,17 @@
status = H5Dclose(dataset);
}
+/****************************************************************************
+ * Spherical harmonics coefficients *
+ ****************************************************************************/
+void h5output_geoid(struct All_variables *E, int cycles)
+{
+
+ /* TODO: write geoid to h5file */
+}
+
+
+
/****************************************************************************
* Create and output /connectivity dataset *
@@ -1914,7 +1929,6 @@
status = set_attribute_float(input, "z_lmantle", E->viscosity.zlm);
status = set_attribute_float(input, "z_cmb", E->viscosity.zcmb);
- status = set_attribute_float(input, "layer_km", E->data.layer_km);
status = set_attribute_float(input, "radius_km", E->data.radius_km);
status = set_attribute_float(input, "scalev", E->data.scalev);
status = set_attribute_float(input, "scalet", E->data.scalet);
@@ -2005,6 +2019,8 @@
status = set_attribute_string(input, "output_format", E->output.format);
status = set_attribute_string(input, "output_optional", E->output.optional);
+ status = set_attribute_int(input, "output_ll_max", E->output.llmax);
+
status = set_attribute_int(input, "verbose", E->control.verbose);
status = set_attribute_int(input, "see_convergence", E->control.print_convergence);
@@ -2070,11 +2086,6 @@
status = set_attribute_double(input, "fi_max", E->control.fi_max);
}
- status = set_attribute_int(input, "ll_max", E->sphere.llmax);
- status = set_attribute_int(input, "nlong", E->sphere.noy);
- status = set_attribute_int(input, "nlati", E->sphere.nox);
- status = set_attribute_int(input, "output_ll_max", E->sphere.output_llmax);
-
/*
* Tracer.inventory
*/
Modified: mc/3D/CitcomS/trunk/lib/Regional_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_obsolete.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Regional_obsolete.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/*
@@ -710,7 +710,7 @@
}
fprintf(fp1,"lmaxx=%.4e lminx=%.4e for %s\n",maxx,minx,filen);
fprintf(fp1," ll mm cos sin \n");
- for (ll=0;ll<=E->sphere.output_llmax;ll++)
+ for (ll=0;ll<=E->output.llmax;ll++)
for(mm=0;mm<=ll;mm++) {
i = E->sphere.hindex[ll][mm];
fprintf(fp1,"%3d %3d %.4e %.4e \n",ll,mm,sphc[i],sphs[i]);
@@ -926,8 +926,8 @@
int m,i,it;
- E->monitor.length_scale = E->data.layer_km/E->mesh.layer[2]; /* km */
- E->monitor.time_scale = pow(E->data.layer_km*1000.0,2.0)/ /* Million years */
+ E->monitor.length_scale = E->data.radius_km/E->mesh.layer[2]; /* km */
+ E->monitor.time_scale = pow(E->data.radius_km*1000.0,2.0)/ /* Million years */
(E->data.therm_diff*3600.0*24.0*365.25*1.0e6);
if ( (ii == 0) || ((ii % E->control.record_every) == 0)
Modified: mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Regional_parallel_related.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -238,6 +238,10 @@
E->lmesh.nsf = E->lmesh.nno/E->lmesh.noz;
E->lmesh.snel = E->lmesh.elx*E->lmesh.ely;
+
+
+
+#if 0
i = cases[E->sphere.caps_per_proc];
E->parallel.nproc_sph[1] = incases3[i].xy[0];
@@ -259,7 +263,10 @@
E->sphere.lexs = E->sphere.lelx * E->parallel.me_loc_sph[1];
E->sphere.leys = E->sphere.lely * E->parallel.me_loc_sph[2];
+#endif
+
+
for(i=E->mesh.levmax;i>=E->mesh.levmin;i--) {
if (E->control.NMULTIGRID||E->control.EMULTIGRID) {
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -114,18 +114,21 @@
}
+#if 0
E->sphere.elx = E->sphere.nox-1;
E->sphere.ely = E->sphere.noy-1;
E->sphere.snel = E->sphere.ely*E->sphere.elx;
E->sphere.nsf = E->sphere.noy*E->sphere.nox;
+#endif
+
/* Scaling from dimensionless units to Millions of years for input velocity
and time, timdir is the direction of time for advection. CPC 6/25/00 */
/* Myr */
- E->data.scalet = (E->data.layer_km*1e3*E->data.layer_km*1e3/E->data.therm_diff)/(1.e6*365.25*24*3600);
+ E->data.scalet = (E->data.radius_km*1e3*E->data.radius_km*1e3/E->data.therm_diff)/(1.e6*365.25*24*3600);
/* cm/yr */
- E->data.scalev = (E->data.layer_km*1e3/E->data.therm_diff)/(100*365.25*24*3600);
+ E->data.scalev = (E->data.radius_km*1e3/E->data.therm_diff)/(100*365.25*24*3600);
E->data.timedir = E->control.Atemp / fabs(E->control.Atemp);
Modified: mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -1,824 +1,339 @@
-/*
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
- *<LicenseText>
- *
- * CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
- * Clint Conrad, Michael Gurnis, and Eun-seo Choi.
- * Copyright (C) 1994-2005, California Institute of Technology.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- *
- *</LicenseText>
- *
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- */
/* Functions relating to the building and use of mesh locations ... */
#include <math.h>
+#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
-#include "parallel_related.h"
+#include <stdlib.h>
-
/* ======================================================================
- ====================================================================== */
+ ====================================================================== */
void set_sphere_harmonics(E)
- struct All_variables *E;
+ struct All_variables *E;
- {
- int node,ll,mm,i,j;
- double multis(),dth,dfi,sqrt_multis();
- void construct_interp_net();
- void compute_sphereh_table();
+{
+ int m,node,ll,mm,i,j;
+ static void compute_sphereh_table();
- E->sphere.sx[1] = (double *) malloc((E->sphere.nsf+1)*sizeof(double));
- E->sphere.sx[2] = (double *) malloc((E->sphere.nsf+1)*sizeof(double));
- E->sphere.int_cap = (int *) malloc((E->sphere.nsf+1)*sizeof(int));
- E->sphere.int_ele = (int *) malloc((E->sphere.nsf+1)*sizeof(int));
-
-/* E->sphere.radius=(float *) malloc((E->sphere.slab_layers+3)*sizeof(float)); */
-
- i=0;
- for (ll=0;ll<=E->sphere.llmax;ll++)
- for (mm=0;mm<=ll;mm++) {
- E->sphere.hindex[ll][mm] = i;
- i++;
+ i=0;
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ E->sphere.hindex[ll][mm] = i;
+ i++;
}
- E->sphere.hindice = i;
- E->sphere.con = (double *)malloc((E->sphere.hindice+3)*sizeof(double));
+ E->sphere.hindice = i;
- for (i=1;i<=E->sphere.lelx;i++)
- E->sphere.tableplm[i]= (double *)malloc((E->sphere.hindice+3)*sizeof(double));
- for (i=1;i<=E->sphere.lely;i++) {
- E->sphere.tablecosf[i]= (double *)malloc((E->sphere.output_llmax+3)*sizeof(double));
- E->sphere.tablesinf[i]= (double *)malloc((E->sphere.output_llmax+3)*sizeof(double));
+ /* spherical harmonic coeff (0=cos, 1=sin)
+ for surface topo, cmb topo and geoid */
+ for (i=0;i<=1;i++) {
+ E->sphere.harm_tpgt[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
+ E->sphere.harm_tpgb[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
+ E->sphere.harm_geoid[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
}
- for (i=1;i<=E->sphere.lnox;i++)
- E->sphere.tableplm_n[i]= (double *)malloc((E->sphere.hindice+3)*sizeof(double));
- for (i=1;i<=E->sphere.lnoy;i++) {
- E->sphere.tablecosf_n[i]= (double *)malloc((E->sphere.output_llmax+3)*sizeof(double));
- E->sphere.tablesinf_n[i]= (double *)malloc((E->sphere.output_llmax+3)*sizeof(double));
- }
- E->sphere.sien = (struct SIEN *) malloc((E->sphere.lsnel+1)*sizeof(struct SIEN));
+ compute_sphereh_table(E);
+ return;
+}
+/* =========================================================
+ expand the field TG into spherical harmonics
+ ========================================================= */
+void sphere_expansion(E,TG,sphc,sphs)
+ struct All_variables *E;
+ float **TG,*sphc,*sphs;
+{
+ int el,nint,d,p,i,m,j,es,mm,ll,rand();
+ //double t,f,sphere_h();
+ void sum_across_surf_sph1();
+ void get_global_1d_shape_fn();
+ struct Shape_function1 M;
+ struct Shape_function1_dA dGamma;
-/* ll = E->sphere.slab_layers+2; */
-
- for (i=0;i<=1;i++) {
- E->sphere.harm_tpgt[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_tpgb[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_velp[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_velt[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_divg[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_vort[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_visc[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
-/* E->sphere.harm_slab[i]=(float*)malloc((ll*E->sphere.hindice+2)*sizeof(float)); */
- }
-
-
-
- for (ll=0;ll<=E->sphere.output_llmax;ll++)
- for (mm=0;mm<=ll;mm++) {
- E->sphere.con[E->sphere.hindex[ll][mm]] =
- sqrt( (2.0-((mm==0)?1.0:0.0))*(2*ll+1)/(4.0*M_PI) )
- *sqrt_multis(ll+mm,ll-mm); /* which is sqrt((ll-mm)!/(ll+mm)!) */
- }
-
- dth = M_PI/E->sphere.elx;
- dfi = 2.0*M_PI/E->sphere.ely;
-
- for (j=1;j<=E->sphere.noy;j++)
- for (i=1;i<=E->sphere.nox;i++) {
- node = i+(j-1)*E->sphere.nox;
- E->sphere.sx[1][node] = dth*(i-1);
- E->sphere.sx[2][node] = dfi*(j-1);
+ for (i=0;i<E->sphere.hindice;i++) {
+ sphc[i] = 0.0;
+ sphs[i] = 0.0;
}
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for (es=1;es<=E->lmesh.snel;es++) {
+ el = es*E->lmesh.elz;
- construct_interp_net(E);
+ get_global_1d_shape_fn(E,el,&M,&dGamma,1,m);
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0; mm<=ll; mm++) {
- compute_sphereh_table(E);
+ p = E->sphere.hindex[ll][mm];
- return;
- }
-/* ======================================================================
- ====================================================================== */
+ for(nint=1;nint<=onedvpoints[E->mesh.nsd];nint++) {
+ for(d=1;d<=onedvpoints[E->mesh.nsd];d++) {
+ j = E->sien[m][es].node[d];
+ sphc[p] += TG[m][E->sien[m][es].node[d]]
+ * E->sphere.tablesplm[m][j][p]
+ * E->sphere.tablescosf[m][j][mm]
+ * E->M.vpt[GMVINDEX(d,nint)]
+ * dGamma.vpt[GMVGAMMA(1,nint)];
+ sphs[p] += TG[m][E->sien[m][es].node[d]]
+ * E->sphere.tablesplm[m][j][p]
+ * E->sphere.tablessinf[m][j][mm]
+ * E->M.vpt[GMVINDEX(d,nint)]
+ * dGamma.vpt[GMVGAMMA(1,nint)];
+ }
+ }
-void sphere_harmonics_layer(E,T,sphc,sphs,iprint,filen)
- struct All_variables *E;
- float **T,*sphc,*sphs;
- int iprint;
- char * filen;
- {
-/* void sphere_expansion(); */
-/* void sphere_interpolate(); */
-/* void print_field_spectral_regular(); */
-/* FILE *fp; */
-/* char output_file[255]; */
-/* int i,node,j,ll,mm,printt,proc_loc; */
-/* float minx,maxx,t,f,rad; */
-/* static int been_here=0; */
-/* float *TG; */
+ } /* end for ll and mm */
-/* rad=180.0/M_PI; */
+ }
-/* maxx=-1.e6; */
-/* minx=1.e6; */
+ sum_across_surf_sph1(E,sphc,sphs);
-/* printt=0; */
+ return;
+}
-/* if(E->parallel.me_loc[3]==E->parallel.nprocz-1 && iprint==1) printt=1; */
-/* if(E->parallel.me_loc[3]==0 && iprint==0) printt=1; */
-/* TG = (float *)malloc((E->sphere.nsf+1)*sizeof(float)); */
+/* ==================================================*/
+/* ==================================================*/
+static void compute_sphereh_table(E)
+ struct All_variables *E;
+{
-/* proc_loc = E->parallel.me_loc[3]; */
+ int m,node,ll,mm,i,j,p;
+ double t,f;
-/* sphere_interpolate(E,T,TG); */
-/* sphere_expansion (E,TG,sphc,sphs); */
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ E->sphere.tablesplm[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
+ E->sphere.tablescosf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
+ E->sphere.tablessinf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
-/* if (printt) */
-/* print_field_spectral_regular(E,TG,sphc,sphs,proc_loc,filen); */
+ for (i=1;i<=E->lmesh.nsf;i++) {
+ E->sphere.tablesplm[m][i]= (double *)malloc((E->sphere.hindice+3)*sizeof(double));
+ E->sphere.tablescosf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
+ E->sphere.tablessinf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
+ }
+ }
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ for (j=1;j<=E->lmesh.nsf;j++) {
+ node = j*E->lmesh.noz;
+ f=E->sx[m][2][node];
+ t=E->sx[m][1][node];
+ for (mm=0;mm<=E->output.llmax;mm++) {
+ E->sphere.tablescosf[m][j][mm] = cos( (double)(mm)*f );
+ E->sphere.tablessinf[m][j][mm] = sin( (double)(mm)*f );
+ }
-/* parallel_process_sync(E); */
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ p = E->sphere.hindex[ll][mm];
+ E->sphere.tablesplm[m][j][p] = modified_plgndr_a(ll,mm,t) ;
+ }
+ }
+ }
-/* free ((void *)TG); */
+ return;
+}
- return;
- }
-
/* ================================================
- compute angle and area
- ================================================*/
+ compute angle and area
+ ================================================*/
void compute_angle_surf_area (E)
- struct All_variables *E;
+ struct All_variables *E;
{
- int es,el,m,i,j,ii,ia[5],lev;
- double aa,y1[4],y2[4],angle[6],xx[4][5],area_sphere_cap();
- void get_angle_sphere_cap();
+ int es,el,m,i,j,ii,ia[5],lev;
+ double aa,y1[4],y2[4],angle[6],xx[4][5],area_sphere_cap();
+ void get_angle_sphere_cap();
+ void parallel_process_termination();
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- ia[1] = 1;
- ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
- ia[3] = E->lmesh.nno-E->lmesh.noz+1;
- ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ ia[1] = 1;
+ ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
+ ia[3] = E->lmesh.nno-E->lmesh.noz+1;
+ ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
- for (i=1;i<=4;i++) {
- xx[1][i] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
- xx[2][i] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
- xx[3][i] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
- }
+ for (i=1;i<=4;i++) {
+ xx[1][i] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
+ xx[2][i] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
+ xx[3][i] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
+ }
- get_angle_sphere_cap(xx,angle);
+ get_angle_sphere_cap(xx,angle);
- for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- E->sphere.angle[m][i] = angle[i];
+ for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ E->sphere.angle[m][i] = angle[i];
- E->sphere.area[m] = area_sphere_cap(angle);
+ E->sphere.area[m] = area_sphere_cap(angle);
- for (lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--)
- for (es=1;es<=E->lmesh.SNEL[lev];es++) {
- el = (es-1)*E->lmesh.ELZ[lev]+1;
- for (i=1;i<=4;i++)
- ia[i] = E->IEN[lev][m][el].node[i];
+ for (lev=E->mesh.levmax;lev>=E->mesh.levmin;lev--)
+ for (es=1;es<=E->lmesh.SNEL[lev];es++) {
+ el = (es-1)*E->lmesh.ELZ[lev]+1;
+ for (i=1;i<=4;i++)
+ ia[i] = E->IEN[lev][m][el].node[i];
- for (i=1;i<=4;i++) {
- xx[1][i] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[2][i] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[3][i] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
- }
+ for (i=1;i<=4;i++) {
+ xx[1][i] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[2][i] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[3][i] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
+ }
- get_angle_sphere_cap(xx,angle);
+ get_angle_sphere_cap(xx,angle);
- for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- E->sphere.angle1[lev][m][i][es] = angle[i];
+ for (i=1;i<=4;i++) /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ E->sphere.angle1[lev][m][i][es] = angle[i];
- E->sphere.area1[lev][m][es] = area_sphere_cap(angle);
+ E->sphere.area1[lev][m][es] = area_sphere_cap(angle);
- /*fprintf(E->fp_out,"lev%d %d %.6e %.6e %.6e %.6e %.6e\n",lev,es,angle[1],angle[2],angle[3],angle[4],E->sphere.area1[lev][m][es]);
-*/
- } /* end for lev and es */
+ fprintf(E->fp_out,"lev%d %d %.6e %.6e %.6e %.6e %.6e\n",lev,es,angle[1],angle[2],angle[3],angle[4],E->sphere.area1[lev][m][es]);
- } /* end for m */
+ } /* end for lev and es */
- return;
- }
+ } /* end for m */
+ return;
+}
+
/* ================================================
- area of spherical rectangle
- ================================================ */
+ area of spherical rectangle
+ ================================================ */
double area_sphere_cap(angle)
- double angle[6];
- {
+ double angle[6];
+{
- double area,a,b,c;
- double area_of_sphere_triag();
+ double area,a,b,c;
+ double area_of_sphere_triag();
- a = angle[1];
- b = angle[2];
- c = angle[5];
- area = area_of_sphere_triag(a,b,c);
+ a = angle[1];
+ b = angle[2];
+ c = angle[5];
+ area = area_of_sphere_triag(a,b,c);
- a = angle[3];
- b = angle[4];
- c = angle[5];
- area += area_of_sphere_triag(a,b,c);
+ a = angle[3];
+ b = angle[4];
+ c = angle[5];
+ area += area_of_sphere_triag(a,b,c);
- return (area);
- }
+ return (area);
+}
/* ================================================
- area of spherical triangle
- ================================================ */
+ area of spherical triangle
+ ================================================ */
double area_of_sphere_triag(a,b,c)
- double a,b,c;
- {
+ double a,b,c;
+{
- double ss,ak,aa,bb,cc,area;
- const double e_16 = 1.0e-16;
- const double two = 2.0;
- const double pt5 = 0.5;
+ double ss,ak,aa,bb,cc,area;
+ const double e_16 = 1.0e-16;
+ const double two = 2.0;
+ const double pt5 = 0.5;
- ss = (a+b+c)*pt5;
- area=0.0;
- a = sin(ss-a);
- b = sin(ss-b);
- c = sin(ss-c);
- ak = a*b*c/sin(ss); /* sin(ss-a)*sin(ss-b)*sin(ss-c)/sin(ss) */
- if(ak<e_16) return (area);
- ak = sqrt(ak);
- aa = two*atan(ak/a);
- bb = two*atan(ak/b);
- cc = two*atan(ak/c);
- area = aa+bb+cc-M_PI;
+ ss = (a+b+c)*pt5;
+ area=0.0;
+ a = sin(ss-a);
+ b = sin(ss-b);
+ c = sin(ss-c);
+ ak = a*b*c/sin(ss); /* sin(ss-a)*sin(ss-b)*sin(ss-c)/sin(ss) */
+ if(ak<e_16) return (area);
+ ak = sqrt(ak);
+ aa = two*atan(ak/a);
+ bb = two*atan(ak/b);
+ cc = two*atan(ak/c);
+ area = aa+bb+cc-M_PI;
- return (area);
- }
+ return (area);
+}
/* =====================================================================
- get the area for given five points (4 nodes for a rectangle and one test node)
- angle [i]: angle bet test node and node i of the rectangle
- angle1[i]: angle bet nodes i and i+1 of the rectangle
- ====================================================================== */
+ get the area for given five points (4 nodes for a rectangle and one test node)
+ angle [i]: angle bet test node and node i of the rectangle
+ angle1[i]: angle bet nodes i and i+1 of the rectangle
+ ====================================================================== */
double area_of_5points(E,lev,m,el,x,ne)
- struct All_variables *E;
- int lev,m,el,ne;
- double x[4];
- {
- int i,es,ia[5];
- double area1,get_angle(),area_of_sphere_triag();
- double xx[4],angle[5],angle1[5];
+ struct All_variables *E;
+ int lev,m,el,ne;
+ double x[4];
+{
+ int i,es,ia[5];
+ double area1,get_angle(),area_of_sphere_triag();
+ double xx[4],angle[5],angle1[5];
for (i=1;i<=4;i++)
- ia[i] = E->IEN[lev][m][el].node[i];
+ ia[i] = E->IEN[lev][m][el].node[i];
es = (el-1)/E->lmesh.ELZ[lev]+1;
for (i=1;i<=4;i++) {
- xx[1] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[2] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
- xx[3] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
- angle[i] = get_angle(x,xx); /* get angle bet (i,j) and other four*/
- angle1[i]= E->sphere.angle1[lev][m][i][es];
- }
+ xx[1] = E->X[lev][m][1][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[2] = E->X[lev][m][2][ia[i]]/E->SX[lev][m][3][ia[1]];
+ xx[3] = E->X[lev][m][3][ia[i]]/E->SX[lev][m][3][ia[1]];
+ angle[i] = get_angle(x,xx); /* get angle bet (i,j) and other four*/
+ angle1[i]= E->sphere.angle1[lev][m][i][es];
+ }
area1 = area_of_sphere_triag(angle[1],angle[2],angle1[1])
- + area_of_sphere_triag(angle[2],angle[3],angle1[2])
- + area_of_sphere_triag(angle[3],angle[4],angle1[3])
- + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
+ + area_of_sphere_triag(angle[2],angle[3],angle1[2])
+ + area_of_sphere_triag(angle[3],angle[4],angle1[3])
+ + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
- return (area1);
- }
+ return (area1);
+}
/* ================================
- get the angle for given four points spherical rectangle
- ================================= */
+ get the angle for given four points spherical rectangle
+ ================================= */
void get_angle_sphere_cap(xx,angle)
- double xx[4][5],angle[6];
- {
+ double xx[4][5],angle[6];
+{
- int i,j,ii;
- double y1[4],y2[4],get_angle();;
+ int i,j,ii;
+ double y1[4],y2[4],get_angle();;
- for (i=1;i<=4;i++) { /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
- for (j=1;j<=3;j++) {
- ii=(i==4)?1:(i+1);
- y1[j] = xx[j][i];
- y2[j] = xx[j][ii];
- }
- angle[i] = get_angle(y1,y2);
- }
+ for (i=1;i<=4;i++) { /* angle1: bet 1 & 2; angle2: bet 2 & 3 ..*/
+ for (j=1;j<=3;j++) {
+ ii=(i==4)?1:(i+1);
+ y1[j] = xx[j][i];
+ y2[j] = xx[j][ii];
+ }
+ angle[i] = get_angle(y1,y2);
+ }
- for (j=1;j<=3;j++) {
- y1[j] = xx[j][1];
- y2[j] = xx[j][3];
- }
+ for (j=1;j<=3;j++) {
+ y1[j] = xx[j][1];
+ y2[j] = xx[j][3];
+ }
- angle[5] = get_angle(y1,y2); /* angle5 for betw 1 and 3: diagonal */
- return;
- }
+ angle[5] = get_angle(y1,y2); /* angle5 for betw 1 and 3: diagonal */
+ return;
+}
/* ================================
- get the angle for given two points
- ================================= */
+ get the angle for given two points
+ ================================= */
double get_angle(x,xx)
- double x[4],xx[4];
- {
- double dist,angle;
- const double pt5 = 0.5;
- const double two = 2.0;
-
- dist=sqrt( (x[1]-xx[1])*(x[1]-xx[1])
- + (x[2]-xx[2])*(x[2]-xx[2])
- + (x[3]-xx[3])*(x[3]-xx[3]) )*pt5;
- angle = asin(dist)*two;
-
- return (angle);
- }
-
-/* ================================================
- for a given node, this routine gives which cap and element
- the node is in.
- ================================================*/
-void construct_interp_net(E)
- struct All_variables *E;
- {
-
-/* int ii,jj,es,i,j,m,el,node; */
-/* int locate_cap(),locate_element(); */
-/* double x[4],t,f; */
-
-/* const int ends=4; */
-
-/* for (i=1;i<=E->sphere.nox;i++) */
-/* for (j=1;j<=E->sphere.noy;j++) { */
-/* node = i+(j-1)*E->sphere.nox; */
-/* E->sphere.int_cap[node]=0; */
-/* E->sphere.int_ele[node]=0; */
-/* } */
-
-
-/* for (i=1;i<=E->sphere.nox;i++) */
-/* for (j=1;j<=E->sphere.noy;j++) { */
-/* node = i+(j-1)*E->sphere.nox */;
-
- /* first find which cap this node (i,j) is in */
-/* t = E->sphere.sx[1][node]; */
-/* f = E->sphere.sx[2][node]; */
-
-/* x[1] = sin(t)*cos(f); */ /* radius does not matter */
-/* x[2] = sin(t)*sin(f); */
-/* x[3] = cos(t); */
-
- /* locate_cap may not work correctly after my change in numbering of caps */
- /* but sphere.int_cap and int_ele are not used anywhere */
-/* m = locate_cap(E,x); */
-/* if (m>0) { */
-/* el = locate_element(E,m,x,node); */ /* bottom element */
-
-/* if (el<=0) { */
-/* fprintf(stderr,"!!! Processor %d cannot find the right element in cap %d\n",E->parallel.me,m); */
-/* parallel_process_termination(); */
-/* } */
-
-/* E->sphere.int_cap[node]=m; */
-/* E->sphere.int_ele[node]=el; */
-
-/* } */
-/* } */ /* end for i and j */
-
-/* parallel_process_sync(E); */
-
- return;
- }
-
-/* ================================================
- locate the cap for node (i,j)
- ================================================*/
-
-int locate_cap(E,x)
- struct All_variables *E;
- double x[4];
- {
-
- int ia[5],i,m,mm;
- double xx[4],angle[5],angle1[5];
- double get_angle();
- double area1,rr;
- const double e_7=1.e-7;
-
- mm = 0;
-
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- ia[1] = 1;
- ia[2] = E->lmesh.noz*E->lmesh.nox-E->lmesh.noz+1;
- ia[3] = E->lmesh.nno-E->lmesh.noz+1;
- ia[4] = ia[3]-E->lmesh.noz*(E->lmesh.nox-1);
-
- for (i=1;i<=4;i++) {
- xx[1] = E->x[m][1][ia[i]]/E->sx[m][3][ia[1]];
- xx[2] = E->x[m][2][ia[i]]/E->sx[m][3][ia[1]];
- xx[3] = E->x[m][3][ia[i]]/E->sx[m][3][ia[1]];
- angle[i] = get_angle(x,xx); /* get angle between (i,j) and other four*/
- angle1[i]=E->sphere.angle[m][i];
- }
-
- area1 = area_of_sphere_triag(angle[1],angle[2],angle1[1])
- + area_of_sphere_triag(angle[2],angle[3],angle1[2])
- + area_of_sphere_triag(angle[3],angle[4],angle1[3])
- + area_of_sphere_triag(angle[4],angle[1],angle1[4]);
-
- if ( fabs ((area1-E->sphere.area[m])/E->sphere.area[m]) <e_7 ) {
- mm = m;
- return (mm);
- }
- }
-
- return (mm);
- }
-
-/* ================================================
- locate the element containing the node (i,j) with coord x.
- The radius is assumed to be 1 in computing the areas.
- NOTE: The returned element el is for the bottom layer.
- ================================================*/
-
-int locate_element(E,m,x,ne)
- struct All_variables *E;
- double x[4];
- int m,ne;
- {
-
- int el_temp,el,es,el_located,level,lev,lev_plus,el_plus,es_plus,i,j,node;
- double area,area1,areamin;
- double area_of_5points();
- const double e_7=1.e-7;
- const double e_6=1.e6;
-
- el_located = 0;
-
-
- level=E->mesh.levmin;
- for (es=1;es<=E->lmesh.SNEL[level];es++) {
-
- el = (es-1)*E->lmesh.ELZ[level]+1;
- area1 = area_of_5points (E,level,m,el,x,ne);
- area = E->sphere.area1[level][m][es];
-
- if(fabs ((area1-area)/area) <e_7 ) {
- for (lev=E->mesh.levmin;lev<E->mesh.levmax;lev++) {
- lev_plus = lev + 1;
- j=1;
- areamin = e_6;
- do {
- el_plus = E->EL[lev][m][el].sub[j];
-
- es_plus = (el_plus-1)/E->lmesh.ELZ[lev_plus]+1;
-
- area1 = area_of_5points(E,lev_plus,m,el_plus,x,ne);
- area = E->sphere.area1[lev_plus][m][es_plus];
-
- if(fabs(area1-area)<areamin) {
- areamin=fabs(area1-area);
- el_temp = el_plus;
- }
- j++;
- } while (j<5 && fabs((area1-area)/area) > e_7);
- el = el_plus;
- /* if exit with ..>e_7, pick the best one*/
- if (fabs((area1-area)/area) > e_7) el = el_temp;
- } /* end for loop lev */
- el_located = el;
- } /* end for if */
-
- if(el_located) break;
- } /* end for es at the coarsest level */
-
- return (el_located);
- }
-
-/* ===============================================================
- interpolate nodal T's within cap m and element el onto node with
- coordinate x[3] which is derived from a regular mesh and within
- the element el. NOTE the radius of x[3] is the inner radius.
- =============================================================== */
-
-float sphere_interpolate_point(E,T,m,el,x,ne)
- struct All_variables *E;
- float **T;
- double x[4];
- int m,el,ne;
- {
-/* double to,fo,y[4],yy[4][5],dxdy[4][4]; */
-/* double a1,b1,c1,d1,a2,b2,c2,d2,a,b,c,xx1,yy1,y1,y2; */
- float ta,t[5];
-/* int es,snode,i,j,node; */
-
-/* const int oned=4; */
-/* const double e_7=1.e-7; */
-/* const double four=4.0; */
-/* const double two=2.0; */
-/* const double one=1.0; */
-/* const double pt25=0.25; */
-
- /* first rotate the coord such that the center of element is
- the pole */
-
-/* es = (el-1)/E->lmesh.elz+1; */
-
-/* to = E->eco[m][el].centre[1]; */
-/* fo = E->eco[m][el].centre[2]; */
-
-/* dxdy[1][1] = cos(to)*cos(fo); */
-/* dxdy[1][2] = cos(to)*sin(fo); */
-/* dxdy[1][3] = -sin(to); */
-/* dxdy[2][1] = -sin(fo); */
-/* dxdy[2][2] = cos(fo); */
-/* dxdy[2][3] = 0.0; */
-/* dxdy[3][1] = sin(to)*cos(fo); */
-/* dxdy[3][2] = sin(to)*sin(fo); */
-/* dxdy[3][3] = cos(to); */
-
-/* for(i=1;i<=oned;i++) { */ /* nodes */
-/* node = E->ien[m][el].node[i]; */
-/* snode = E->sien[m][es].node[i]; */
-/* t[i] = T[m][snode]; */
-/* for (j=1;j<=E->mesh.nsd;j++) */
-/* yy[j][i] = E->x[m][1][node]*dxdy[j][1] */
-/* + E->x[m][2][node]*dxdy[j][2] */
-/* + E->x[m][3][node]*dxdy[j][3]; */
-/* } */
-
-/* for (j=1;j<=E->mesh.nsd;j++) */
-/* y[j] = x[1]*dxdy[j][1] + x[2]*dxdy[j][2] + x[3]*dxdy[j][3]; */
-
- /* then for node y, determine its coordinates xx1,yy1
- in the parental element in the isoparametric element system*/
-
-/* a1 = yy[1][1] + yy[1][2] + yy[1][3] + yy[1][4]; */
-/* b1 = yy[1][3] + yy[1][2] - yy[1][1] - yy[1][4]; */
-/* c1 = yy[1][3] + yy[1][1] - yy[1][2] - yy[1][4]; */
-/* d1 = yy[1][3] + yy[1][4] - yy[1][1] - yy[1][2]; */
-/* a2 = yy[2][1] + yy[2][2] + yy[2][3] + yy[2][4]; */
-/* b2 = yy[2][3] + yy[2][2] - yy[2][1] - yy[2][4]; */
-/* c2 = yy[2][3] + yy[2][1] - yy[2][2] - yy[2][4]; */
-/* d2 = yy[2][3] + yy[2][4] - yy[2][1] - yy[2][2]; */
-
-/* a = d2*c1; */
-/* b = a2*c1+b1*d2-d1*c2-d1*b2-four*c1*y[2]; */
-/* c=four*c2*y[1]-c2*a1-a1*b2+four*b2*y[1]-four*b1*y[2]+a2*b1; */
-
-/* if (fabs(a)<e_7) { */
-/* yy1 = -c/b; */
-/* xx1 = (four*y[1]-a1-d1*yy1)/(b1+c1*yy1); */
-/* } */
-/* else { */
-/* y1= (-b+sqrt(b*b-four*a*c))/(two*a); */
-/* y2= (-b-sqrt(b*b-four*a*c))/(two*a); */
-/* if (fabs(y1)>fabs(y2)) */
-/* yy1 = y2; */
-/* else */
-/* yy1 = y1; */
-/* xx1 = (four*y[1]-a1-d1*yy1)/(b1+c1*yy1); */
-/* } */
-
- /* now we can calculate T at x[4] using shape function */
-
-/* ta = ((one-xx1)*(one-yy1)*t[1]+(one+xx1)*(one-yy1)*t[2]+ */
-/* (one+xx1)*(one+yy1)*t[3]+(one-xx1)*(one+yy1)*t[4])*pt25; */
-
-/*if(fabs(xx1)>1.5 || fabs(yy1)>1.5)fprintf(E->fp_out,"ME= %d %d %d %g %g %g %g %g %g %g\n",ne,m,es,t[1],t[2],t[3],t[4],ta,xx1,yy1);
- */
- return (ta);
- }
-
-/* ===================================================================
- do the interpolation on sphere for data T, which is needed for both
- spherical harmonic expansion and graphics
- =================================================================== */
-
-void sphere_interpolate(E,T,TG)
- struct All_variables *E;
- float **T,*TG;
- {
-
-/* float sphere_interpolate_point(); */
-/* void gather_TG_to_me0(); */
-
-/* int ii,jj,es,i,j,m,el,node; */
-/* double x[4],t,f; */
-
-/* const int ends=4; */
-
-/* for (i=1;i<=E->sphere.nox;i++) */
-/* for (j=1;j<=E->sphere.noy;j++) { */
-/* node = i+(j-1)*E->sphere.nox; */
-/* TG[node] = 0.0; */
- /* first find which cap this node (i,j) is in */
-
-/* m = E->sphere.int_cap[node]; */
-/* el = E->sphere.int_ele[node]; */
-
-/* if (m>0 && el>0) { */
-/* t = E->sphere.sx[1][node]; */
-/* f = E->sphere.sx[2][node]; */
-
-/* x[1] = E->sx[1][3][1]*sin(t)*cos(f); */
-/* x[2] = E->sx[1][3][1]*sin(t)*sin(f); */
-/* x[3] = E->sx[1][3][1]*cos(t); */
-
-/* TG[node] = sphere_interpolate_point(E,T,m,el,x,node); */
-
-/* } */
-
-/* } */
-
-/* gather_TG_to_me0(E,TG); */
-
- return;
- }
-
-/* =========================================================
- ========================================================= */
-void sphere_expansion(E,TG,sphc,sphs)
- struct All_variables *E;
- float *TG,*sphc,*sphs;
- {
-/* int p,i,j,es,mm,ll,rand(); */
-/* double temp,area,t,f,sphere_h(); */
-/* const double pt25=0.25; */
-/* static int been_here=0; */
-/* void sum_across_surf_sph1(); */
-/* void sum_across_surf_sph(); */
-
-/* for (i=0;i<E->sphere.hindice;i++) { */
-/* sphc[i] = 0.0; */
-/* sphs[i] = 0.0; */
-/* } */
-
-/* area = 2.0*M_PI*M_PI/(E->sphere.elx*E->sphere.ely); */
-
-/* for (ll=0;ll<=E->sphere.output_llmax;ll++) */
-/* for (mm=0; mm<=ll; mm++) { */
-
-/* p = E->sphere.hindex[ll][mm]; */
-
-/* for (j=1;j<=E->sphere.lely;j++) */
-/* for (i=1;i<=E->sphere.lelx;i++) { */
-
-/* es = i+(j-1)*E->sphere.lelx; */
-
-/* temp=pt25*(TG[E->sphere.sien[es].node[1]] */
-/* +TG[E->sphere.sien[es].node[2]] */
-/* +TG[E->sphere.sien[es].node[3]] */
-/* +TG[E->sphere.sien[es].node[4]]); */
-
-/* sphc[p]+=temp*E->sphere.tableplm[i][p]*E->sphere.tablecosf[j][mm]*E->sphere.tablesint[i]; */
-/* sphs[p]+=temp*E->sphere.tableplm[i][p]*E->sphere.tablesinf[j][mm]*E->sphere.tablesint[i]; */
-
-/* } */
-
-/* sphc[p] *= area; */
-/* sphs[p] *= area; */
-
-/* } */ /* end for ll and mm */
-
-/* sum_across_surf_sph1(E,sphc,sphs); */
-
- return;
- }
-
- /* =========================================================== */
-void inv_sphere_harmonics(E,sphc,sphs,TG,proc_loc)
- struct All_variables *E;
- float *TG,*sphc,*sphs;
- int proc_loc;
- {
-/* int k,ll,mm,node,i,j,p,noz,snode; */
-/* float t1,f1,rad; */
-/* void gather_TG_to_me0(); */
-
-/* if (E->parallel.me_loc[3]==proc_loc) { */
-
-/* for (j=1;j<=E->sphere.noy;j++) */
-/* for (i=1;i<=E->sphere.nox;i++) { */
-/* node = i + (j-1)*E->sphere.nox; */
-/* TG[node]=0.0; */
-/* } */
-
-/* for (ll=0;ll<=E->sphere.output_llmax;ll++) */
-/* for (mm=0;mm<=ll;mm++) { */
-
-/* p = E->sphere.hindex[ll][mm]; */
-
-/* for (i=1;i<=E->sphere.lnox;i++) */
-/* for (j=1;j<=E->sphere.lnoy;j++) { */
-/* node = i + E->sphere.lexs + (j+E->sphere.leys-1)*E->sphere.nox; */
-/* TG[node] += */
-/* (sphc[p]*E->sphere.tableplm_n[i][p]*E->sphere.tablecosf_n[j][mm] */
-/* +sphs[p]*E->sphere.tableplm_n[i][p]*E->sphere.tablesinf_n[j][mm]); */
-/* } */
-/* } */
-
-/* gather_TG_to_me0(E,TG); */
-
-/* } */
-
-/* parallel_process_sync(E); */
-
- return;
- }
-
-/* ==================================================*/
-void compute_sphereh_table(E)
-struct All_variables *E;
+ double x[4],xx[4];
{
+ double dist,angle;
+ const double pt5 = 0.5;
+ const double two = 2.0;
-int rr,node,ends,ll,mm,es,i,j,p;
-double t,f,plgndr_a(),modified_plgndr_a();
-const double pt25=0.25;
+ dist=sqrt( (x[1]-xx[1])*(x[1]-xx[1])
+ + (x[2]-xx[2])*(x[2]-xx[2])
+ + (x[3]-xx[3])*(x[3]-xx[3]) )*pt5;
+ angle = asin(dist)*two;
- ends = 4;
+ return (angle);
+}
-/* for (j=1;j<=E->sphere.lely;j++) */
-/* for (i=1;i<=E->sphere.lelx;i++) { */
-/* es = i+(j-1)*E->sphere.lelx; */
-/* node = E->sphere.lexs + i + (E->sphere.leys+j-1)*E->sphere.nox; */
-/* for (rr=1;rr<=ends;rr++) */
-/* E->sphere.sien[es].node[rr] = node */
-/* + offset[rr].vector[1] */
-/* + offset[rr].vector[2]*E->sphere.nox; */
-/* } */
-/* for (j=1;j<=E->sphere.lely;j++) { */
-/* es = 1+(j-1)*E->sphere.lelx; */
-/* f=pt25*(E->sphere.sx[2][E->sphere.sien[es].node[1]] */
-/* +E->sphere.sx[2][E->sphere.sien[es].node[2]] */
-/* +E->sphere.sx[2][E->sphere.sien[es].node[3]] */
-/* +E->sphere.sx[2][E->sphere.sien[es].node[4]]); */
-/* for (mm=0;mm<=E->sphere.output_llmax;mm++) { */
-/* E->sphere.tablecosf[j][mm] = cos( (double)(mm)*f ); */
-/* E->sphere.tablesinf[j][mm] = sin( (double)(mm)*f ); */
-/* } */
-/* } */
+#if 0
-/* for (i=1;i<=E->sphere.lelx;i++) { */
-/* es = i+(1-1)*E->sphere.lelx; */
-/* t=pt25*(E->sphere.sx[1][E->sphere.sien[es].node[1]] */
-/* +E->sphere.sx[1][E->sphere.sien[es].node[2]] */
-/* +E->sphere.sx[1][E->sphere.sien[es].node[3]] */
-/* +E->sphere.sx[1][E->sphere.sien[es].node[4]]); */
-/* E->sphere.tablesint[i] = sin(t); */
-/* for (ll=0;ll<=E->sphere.output_llmax;ll++) */
-/* for (mm=0;mm<=ll;mm++) { */
-/* p = E->sphere.hindex[ll][mm]; */
-/* E->sphere.tableplm[i][p] = modified_plgndr_a(ll,mm,t) ; */
-/* } */
-/* } */
-/* for (j=1;j<=E->sphere.lnoy;j++) { */
-/* node = E->sphere.lexs + 1 + (E->sphere.leys+j-1)*E->sphere.nox; */
-/* f=E->sphere.sx[2][node]; */
-/* for (mm=0;mm<=E->sphere.output_llmax;mm++) { */
-/* E->sphere.tablecosf_n[j][mm] = cos( (double)(mm)*f ); */
-/* E->sphere.tablesinf_n[j][mm] = sin( (double)(mm)*f ); */
-/* } */
-/* } */
-
-/* for (i=1;i<=E->sphere.lnox;i++) { */
-/* node = E->sphere.lexs + i + (E->sphere.leys+1-1)*E->sphere.nox; */
-/* t=E->sphere.sx[1][node]; */
-/* for (ll=0;ll<=E->sphere.output_llmax;ll++) */
-/* for (mm=0;mm<=ll;mm++) { */
-/* p = E->sphere.hindex[ll][mm]; */
-/* E->sphere.tableplm_n[i][p] = modified_plgndr_a(ll,mm,t) ; */
-/* } */
-/* } */
-
- return;
- }
+#endif
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2006-11-16 01:52:25 UTC (rev 5296)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <stdio.h>
@@ -77,23 +77,23 @@
void get_STD_freesurf(struct All_variables *E,float **freesurf)
{
- int node,snode,m;
+ int node,snode,m;
- if (E->parallel.me_loc[3]==E->parallel.nprocz-1)
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(snode=1;snode<=E->lmesh.nsf;snode++) {
- node = E->surf_node[m][snode];
- /*freesurf[m][snode] += 0.5*(E->sphere.cap[m].V[3][node]+E->sphere.cap[m].Vprev[3][node])*E->advection.timestep;*/
- freesurf[m][snode] += E->sphere.cap[m].V[3][node]*E->advection.timestep;
- }
- return;
+ if (E->parallel.me_loc[3]==E->parallel.nprocz-1)
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(snode=1;snode<=E->lmesh.nsf;snode++) {
+ node = E->surf_node[m][snode];
+ /*freesurf[m][snode] += 0.5*(E->sphere.cap[m].V[3][node]+E->sphere.cap[m].Vprev[3][node])*E->advection.timestep;*/
+ freesurf[m][snode] += E->sphere.cap[m].V[3][node]*E->advection.timestep;
+ }
+ return;
}
void allocate_STD_mem(struct All_variables *E,
- float** SXX, float** SYY, float** SZZ,
- float** SXY, float** SXZ, float** SZY,
- float** divv, float** vorv)
+ float** SXX, float** SYY, float** SZZ,
+ float** SXY, float** SXZ, float** SZY,
+ float** divv, float** vorv)
{
int m, i;
@@ -125,9 +125,9 @@
void free_STD_mem(struct All_variables *E,
- float** SXX, float** SYY, float** SZZ,
- float** SXY, float** SXZ, float** SZY,
- float** divv, float** vorv)
+ float** SXX, float** SYY, float** SZZ,
+ float** SXY, float** SXZ, float** SZY,
+ float** divv, float** vorv)
{
int m;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -177,9 +177,9 @@
void compute_nodal_stress(struct All_variables *E,
- float** SXX, float** SYY, float** SZZ,
- float** SXY, float** SXZ, float** SZY,
- float** divv, float** vorv)
+ float** SXX, float** SYY, float** SZZ,
+ float** SXY, float** SXZ, float** SZY,
+ float** divv, float** vorv)
{
void get_global_shape_fn();
void velo_from_element();
@@ -216,57 +216,57 @@
velo_from_element(E,VV,m,e,sphere_key);
for(j=1;j<=vpts;j++) {
- pre[j] = E->EVi[m][(e-1)*vpts+j]*dOmega.vpt[j];
- Vxyz[1][j] = 0.0;
- Vxyz[2][j] = 0.0;
- Vxyz[3][j] = 0.0;
- Vxyz[4][j] = 0.0;
- Vxyz[5][j] = 0.0;
- Vxyz[6][j] = 0.0;
- Vxyz[7][j] = 0.0;
- Vxyz[8][j] = 0.0;
+ pre[j] = E->EVi[m][(e-1)*vpts+j]*dOmega.vpt[j];
+ Vxyz[1][j] = 0.0;
+ Vxyz[2][j] = 0.0;
+ Vxyz[3][j] = 0.0;
+ Vxyz[4][j] = 0.0;
+ Vxyz[5][j] = 0.0;
+ Vxyz[6][j] = 0.0;
+ Vxyz[7][j] = 0.0;
+ Vxyz[8][j] = 0.0;
}
for(i=1;i<=ends;i++) {
- tww[i] = 0.0;
- for(j=1;j<=vpts;j++)
- tww[i] += dOmega.vpt[j] * g_point[j].weight[E->mesh.nsd-1]
- * E->N.vpt[GNVINDEX(i,j)];
+ tww[i] = 0.0;
+ for(j=1;j<=vpts;j++)
+ tww[i] += dOmega.vpt[j] * g_point[j].weight[E->mesh.nsd-1]
+ * E->N.vpt[GNVINDEX(i,j)];
}
for(j=1;j<=vpts;j++) {
- for(i=1;i<=ends;i++) {
- Vxyz[1][j]+=( VV[1][i]*GNx.vpt[GNVXINDEX(0,i,j)]
- + VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
- Vxyz[2][j]+=( (VV[2][i]*GNx.vpt[GNVXINDEX(1,i,j)]
- + VV[1][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j]))/sin(rtf[1][j])
- + VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
- Vxyz[3][j]+= VV[3][i]*GNx.vpt[GNVXINDEX(2,i,j)];
+ for(i=1;i<=ends;i++) {
+ Vxyz[1][j]+=( VV[1][i]*GNx.vpt[GNVXINDEX(0,i,j)]
+ + VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
+ Vxyz[2][j]+=( (VV[2][i]*GNx.vpt[GNVXINDEX(1,i,j)]
+ + VV[1][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j]))/sin(rtf[1][j])
+ + VV[3][i]*E->N.vpt[GNVINDEX(i,j)] )*rtf[3][j];
+ Vxyz[3][j]+= VV[3][i]*GNx.vpt[GNVXINDEX(2,i,j)];
- Vxyz[4][j]+=( (VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)]
- - VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j]))/sin(rtf[1][j])
- + VV[2][i]*GNx.vpt[GNVXINDEX(0,i,j)])*rtf[3][j];
- Vxyz[5][j]+=VV[1][i]*GNx.vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
- *GNx.vpt[GNVXINDEX(0,i,j)]-VV[1][i]*E->N.vpt[GNVINDEX(i,j)]);
- Vxyz[6][j]+=VV[2][i]*GNx.vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
- *GNx.vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j])-VV[2][i]*E->N.vpt[GNVINDEX(i,j)]);
- Vxyz[7][j]+=rtf[3][j] * (
- VV[1][i]*GNx.vpt[GNVXINDEX(0,i,j)]
- + VV[1][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])/sin(rtf[1][j])
- + VV[2][i]*GNx.vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j]) );
- Vxyz[8][j]+=rtf[3][j]/sin(rtf[1][j])*
- ( VV[2][i]*GNx.vpt[GNVXINDEX(0,i,j)]*sin(rtf[1][j])
- + VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])
- - VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)] );
- }
- Sxx += 2.0 * pre[j] * Vxyz[1][j];
- Syy += 2.0 * pre[j] * Vxyz[2][j];
- Szz += 2.0 * pre[j] * Vxyz[3][j];
- Sxy += pre[j] * Vxyz[4][j];
- Sxz += pre[j] * Vxyz[5][j];
- Szy += pre[j] * Vxyz[6][j];
- div += Vxyz[7][j]*dOmega.vpt[j];
- vor += Vxyz[8][j]*dOmega.vpt[j];
+ Vxyz[4][j]+=( (VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)]
+ - VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j]))/sin(rtf[1][j])
+ + VV[2][i]*GNx.vpt[GNVXINDEX(0,i,j)])*rtf[3][j];
+ Vxyz[5][j]+=VV[1][i]*GNx.vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
+ *GNx.vpt[GNVXINDEX(0,i,j)]-VV[1][i]*E->N.vpt[GNVINDEX(i,j)]);
+ Vxyz[6][j]+=VV[2][i]*GNx.vpt[GNVXINDEX(2,i,j)] + rtf[3][j]*(VV[3][i]
+ *GNx.vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j])-VV[2][i]*E->N.vpt[GNVINDEX(i,j)]);
+ Vxyz[7][j]+=rtf[3][j] * (
+ VV[1][i]*GNx.vpt[GNVXINDEX(0,i,j)]
+ + VV[1][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])/sin(rtf[1][j])
+ + VV[2][i]*GNx.vpt[GNVXINDEX(1,i,j)]/sin(rtf[1][j]) );
+ Vxyz[8][j]+=rtf[3][j]/sin(rtf[1][j])*
+ ( VV[2][i]*GNx.vpt[GNVXINDEX(0,i,j)]*sin(rtf[1][j])
+ + VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])
+ - VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)] );
+ }
+ Sxx += 2.0 * pre[j] * Vxyz[1][j];
+ Syy += 2.0 * pre[j] * Vxyz[2][j];
+ Szz += 2.0 * pre[j] * Vxyz[3][j];
+ Sxy += pre[j] * Vxyz[4][j];
+ Sxz += pre[j] * Vxyz[5][j];
+ Szy += pre[j] * Vxyz[6][j];
+ div += Vxyz[7][j]*dOmega.vpt[j];
+ vor += Vxyz[8][j]*dOmega.vpt[j];
}
Sxx /= E->eco[m][e].area;
@@ -283,15 +283,15 @@
Syy -= E->P[m][e]; /* add the pressure term */
for(i=1;i<=ends;i++) {
- node = E->ien[m][e].node[i];
- SZZ[m][node] += tww[i] * Szz;
- SXX[m][node] += tww[i] * Sxx;
- SYY[m][node] += tww[i] * Syy;
- SXY[m][node] += tww[i] * Sxy;
- SXZ[m][node] += tww[i] * Sxz;
- SZY[m][node] += tww[i] * Szy;
- divv[m][node]+= tww[i] * div;
- vorv[m][node]+= tww[i] * vor;
+ node = E->ien[m][e].node[i];
+ SZZ[m][node] += tww[i] * Szz;
+ SXX[m][node] += tww[i] * Sxx;
+ SYY[m][node] += tww[i] * Syy;
+ SXY[m][node] += tww[i] * Sxy;
+ SXZ[m][node] += tww[i] * Sxz;
+ SZY[m][node] += tww[i] * Szy;
+ divv[m][node]+= tww[i] * div;
+ vorv[m][node]+= tww[i] * vor;
}
} /* end for el */
@@ -350,58 +350,273 @@
int m, i, j, k, n, d;
const unsigned sbc_flag[4] = {0, SBX, SBY, SBZ};
const int stress_index[4][4] = { {0, 0, 0, 0},
- {0, 1, 4, 5},
- {0, 4, 2, 6},
- {0, 5, 6, 3} };
+ {0, 1, 4, 5},
+ {0, 4, 2, 6},
+ {0, 5, 6, 3} };
if(E->control.side_sbcs) {
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(i=1; i<=E->lmesh.noy; i++)
for(j=1; j<=E->lmesh.nox; j++)
- for(k=1; k<=E->lmesh.noz; k++) {
- n = k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
- for(d=1; d<=E->mesh.nsd; d++)
- if(E->node[m][n] & sbc_flag[d]) {
- if(i==1)
- E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sbc.SB[m][SIDE_WEST][d][ E->sbc.node[m][n] ];
- if(i==E->lmesh.noy)
- E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sbc.SB[m][SIDE_EAST][d][ E->sbc.node[m][n] ];
- if(j==1)
- E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sbc.SB[m][SIDE_NORTH][d][ E->sbc.node[m][n] ];
- if(j==E->lmesh.nox)
- E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sbc.SB[m][SIDE_SOUTH][d][ E->sbc.node[m][n] ];
- if(k==1)
- E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sbc.SB[m][SIDE_BOTTOM][d][ E->sbc.node[m][n] ];
- if(k==E->lmesh.noz)
- E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sbc.SB[m][SIDE_TOP][d][ E->sbc.node[m][n] ];
- }
- }
+ for(k=1; k<=E->lmesh.noz; k++) {
+ n = k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
+ for(d=1; d<=E->mesh.nsd; d++)
+ if(E->node[m][n] & sbc_flag[d]) {
+ if(i==1)
+ E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sbc.SB[m][SIDE_WEST][d][ E->sbc.node[m][n] ];
+ if(i==E->lmesh.noy)
+ E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sbc.SB[m][SIDE_EAST][d][ E->sbc.node[m][n] ];
+ if(j==1)
+ E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sbc.SB[m][SIDE_NORTH][d][ E->sbc.node[m][n] ];
+ if(j==E->lmesh.nox)
+ E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sbc.SB[m][SIDE_SOUTH][d][ E->sbc.node[m][n] ];
+ if(k==1)
+ E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sbc.SB[m][SIDE_BOTTOM][d][ E->sbc.node[m][n] ];
+ if(k==E->lmesh.noz)
+ E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sbc.SB[m][SIDE_TOP][d][ E->sbc.node[m][n] ];
+ }
+ }
} else {
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(i=1; i<=E->lmesh.noy; i++)
for(j=1; j<=E->lmesh.nox; j++)
- for(k=1; k<=E->lmesh.noz; k++) {
- n = k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
- for(d=1; d<=E->mesh.nsd; d++)
- if(E->node[m][n] & sbc_flag[d]) {
- if(i==1 || i==E->lmesh.noy)
- E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sphere.cap[m].VB[d][n];
- if(j==1 || j==E->lmesh.nox)
- E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sphere.cap[m].VB[d][n];
- if(k==1 || k==E->lmesh.noz)
- E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sphere.cap[m].VB[d][n];
- }
- }
+ for(k=1; k<=E->lmesh.noz; k++) {
+ n = k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
+ for(d=1; d<=E->mesh.nsd; d++)
+ if(E->node[m][n] & sbc_flag[d]) {
+ if(i==1 || i==E->lmesh.noy)
+ E->gstress[m][(n-1)*6+stress_index[d][2]] = E->sphere.cap[m].VB[d][n];
+ if(j==1 || j==E->lmesh.nox)
+ E->gstress[m][(n-1)*6+stress_index[d][1]] = E->sphere.cap[m].VB[d][n];
+ if(k==1 || k==E->lmesh.noz)
+ E->gstress[m][(n-1)*6+stress_index[d][3]] = E->sphere.cap[m].VB[d][n];
+ }
+ }
}
}
/* ===================================================================
=================================================================== */
+void compute_geoid(E, harm_geoid, harm_tpgt, harm_tpgb)
+ struct All_variables *E;
+ float *harm_geoid[2], *harm_tpgt[2], *harm_tpgb[2];
+{
+ static void geoid_from_buoyancy();
+ static void geoid_from_topography();
+ int ll, mm, p;
+
+ for (ll=1;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ p = E->sphere.hindex[ll][mm];
+ harm_geoid[0][p] = 0.0;
+ harm_geoid[1][p] = 0.0;
+ }
+
+ geoid_from_buoyancy(E, harm_geoid);
+ geoid_from_topography(E, harm_geoid, harm_tpgt, harm_tpgb);
+}
+
+
+static void geoid_from_buoyancy(E, harm_geoid)
+ struct All_variables *E;
+ float *harm_geoid[2];
+{
+ int m,k,ll,mm,node,i,j,p,noz,snode;
+ float *TT[NCS],radius,*geoid[2],dlayer,con1,con2,scaling2,scaling;
+ void sphere_expansion();
+ void sum_across_depth_sph1();
+
+ /* scale for buoyancy */
+ scaling2 = -E->data.therm_exp*E->data.ref_temperature*E->data.density/E->control.Atemp;
+ /* scale for geoid */
+ scaling = 1.0e6*4.0*M_PI*E->data.radius_km*E->data.radius_km*
+ E->data.grav_const/E->data.grav_acc;
+
+ /* buoyancy of one layer */
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ TT[m] = (float *) malloc ((E->lmesh.nsf+1)*sizeof(float));
+
+ /* sin coeff */
+ geoid[0] = (float*)malloc((E->sphere.hindice+1)*sizeof(float));
+ /* cos coeff */
+ geoid[1] = (float*)malloc((E->sphere.hindice+1)*sizeof(float));
+
+ /* loop over each layer */
+ for(k=1;k<E->lmesh.noz;k++) {
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.noy;i++)
+ for(j=1;j<=E->lmesh.nox;j++) {
+ node=k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
+ p = j + (i-1)*E->lmesh.nox;
+ TT[m][p] = (E->buoyancy[m][node]+E->buoyancy[m][node+1])*0.5*scaling2;
+ }
+
+ /* expand TT into spherical harmonics */
+ sphere_expansion(E,TT,geoid[0],geoid[1]);
+
+ /* thickness of the layer */
+ dlayer = (E->sx[1][3][k+1]-E->sx[1][3][k]);
+
+ /* radius of the layer */
+ radius = (E->sx[1][3][k+1]+E->sx[1][3][k])*0.5;
+
+ /* contribution of buoyancy at this layer */
+ for (ll=1;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ con1 = scaling/(2.0*ll+1.0)*dlayer;
+ con2 = pow(radius,((double)(ll+2)));
+
+ p = E->sphere.hindex[ll][mm];
+ harm_geoid[0][p] += con1*con2*geoid[0][p];
+ harm_geoid[1][p] += con1*con2*geoid[1][p];
+ }
+
+ //if(E->parallel.me==0) fprintf(stderr,"layer %d %.5e %g %g %g\n",k,radius,dlayer,con1,con2);
+ }
+
+ /* accumulate geoid from all layers to the surface (top processors) */
+ sum_across_depth_sph1(E, harm_geoid[0], harm_geoid[1]);
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ free ((void *)TT[m]);
+
+ free ((void *)geoid[0]);
+ free ((void *)geoid[1]);
+ return;
+}
+
+
+
+static void geoid_from_topography(E, harm_geoid, harm_tpgt, harm_tpgb)
+ struct All_variables *E;
+ float *harm_geoid[2], *harm_tpgt[2], *harm_tpgb[2];
+{
+
+ float *geoid[2],con1,con2,scaling,den_contrast1,den_contrast2,stress_scaling,topo_scaling1,topo_scaling2;
+ int i,j,k,ll,mm,s;
+ void sum_across_depth_sph1();
+ void sphere_expansion();
+
+ geoid[0] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
+ geoid[1] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
+
+ for (i=0;i<E->sphere.hindice;i++) {
+ geoid[0][i] = 0.0;
+ geoid[1][i] = 0.0;
+ }
+
+ stress_scaling = E->data.ref_viscosity*E->data.therm_diff/
+ (E->data.radius_km*E->data.radius_km*1e6);
+
+ /* density across surface */
+ den_contrast1 = (E->data.density-E->data.density_above);
+ /* density across CMB */
+ den_contrast2 = (E->data.density_below-E->data.density);
+
+ /* scale for surface and CMB topo */
+ topo_scaling1 = stress_scaling / (den_contrast1*E->data.grav_acc);
+ topo_scaling2 = stress_scaling / (den_contrast2*E->data.grav_acc);
+
+ /* scale for geoid */
+ scaling = 1.0e3*4.0*M_PI*E->data.radius_km*
+ E->data.grav_const/E->data.grav_acc;
+
+
+ if (E->parallel.me_loc[3]==E->parallel.nprocz-1) {
+ /* dimensionalize surface topography and
+ expand into spherical harmonics */
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ for(i=1;i<=E->lmesh.nsf;i++) {
+ E->slice.tpg[j][i] = E->slice.tpg[j][i]*topo_scaling1;
+ }
+ sphere_expansion(E,E->slice.tpg,E->sphere.harm_tpgt[0],
+ E->sphere.harm_tpgt[1]);
+ }
+ //print_field_spectral_regular(E,TL,E->sphere.harm_tpgt[0],E->sphere.harm_tpgt[1],E->parallel.nprocz-1,0,"tpgt");
+
+
+
+
+ if (E->parallel.me_loc[3]==0) {
+ /* dimensionalize CMB topography and
+ expand into spherical harmonics */
+ for(j=1;j<=E->sphere.caps_per_proc;j++)
+ for(i=1;i<=E->lmesh.nsf;i++) {
+ E->slice.tpgb[j][i] = E->slice.tpgb[j][i]*topo_scaling2;
+ }
+ sphere_expansion(E,E->slice.tpgb,E->sphere.harm_tpgb[0],
+ E->sphere.harm_tpgb[1]);
+ }
+ //print_field_spectral_regular(E,TL,E->sphere.harm_tpgb[0],E->sphere.harm_tpgb[1],0,0,"tpgb");
+
+
+
+
+
+ if (E->parallel.me_loc[3] == E->parallel.nprocz-1)
+ for (ll=2;ll<=E->output.llmax;ll++) {
+ con1 = scaling/(2.0*ll+1.0);
+ for (mm=0;mm<=ll;mm++) {
+ i = E->sphere.hindex[ll][mm];
+ geoid[0][i]+= E->sphere.harm_tpgt[0][i]*con1*den_contrast1;
+ geoid[1][i]+= E->sphere.harm_tpgt[1][i]*con1*den_contrast1;
+ }
+ }
+
+
+
+
+ if (E->parallel.me_loc[3] == 0)
+ for (ll=2;ll<=E->output.llmax;ll++) {
+ con1 = scaling/(2.0*ll+1.0);
+ con2 = pow(E->sphere.ri,((double)(ll+2)));
+ for (mm=0;mm<=ll;mm++) {
+ i = E->sphere.hindex[ll][mm];
+ geoid[0][i]+= E->sphere.harm_tpgb[0][i]*con1*con2*den_contrast2;
+ geoid[1][i]+= E->sphere.harm_tpgb[1][i]*con1*con2*den_contrast2;
+ }
+ }
+
+ /* accumulate geoid to the surface (top processors) */
+ sum_across_depth_sph1(E,geoid[0],geoid[1]);
+
+ // add to the geoid coeff
+ for (ll=2;ll<=E->output.llmax;ll++)
+ for (mm=0;mm<=ll;mm++) {
+ i = E->sphere.hindex[ll][mm];
+ harm_geoid[0][i]+= geoid[0][i];
+ harm_geoid[1][i]+= geoid[1][i];
+ }
+
+
+
+#if 0
+ E->sphere.harm_geoid[0][0]=0.0;
+ E->sphere.harm_geoid[1][0]=0.0;
+ E->sphere.harm_geoid[0][1]=0.0;
+ E->sphere.harm_geoid[1][1]=0.0;
+ E->sphere.harm_geoid[0][2]=0.0;
+ E->sphere.harm_geoid[1][2]=0.0; /* zero the 00 10 11 */
+
+ inv_sphere_harmonics(E,E->sphere.harm_geoid[0],E->sphere.harm_geoid[1],TL,E->parallel.nprocz-1);
+
+ print_field_spectral_regular(E,TL,E->sphere.harm_geoid[0],E->sphere.harm_geoid[1],E->parallel.nprocz-1,1,"geoid");
+#endif
+
+
+
+ free ((void *)geoid[0]);
+ free ((void *)geoid[1]);
+
+ return;
+}
+
/* ===================================================================
Consistent boundary flux method for stress ... Zhong,Gurnis,Hulbert
@@ -508,23 +723,23 @@
/* resb[m] -= eltkb[ends*dims*m+l] * eub[l]; */
/* } */
- /* Put relevant (vertical & surface) parts of element residual into surface residual */
+ /* Put relevant (vertical & surface) parts of element residual into surface residual */
/* for(m=1;m<=ends;m++) { */ /* for bottom elements */
/* switch (m) { */
/* case 2: */
-/* RL[j][E->sien[j][els].node[1]] += resb[(m-1)*dims+1]; */
-/* break; */
+/* RL[j][E->sien[j][els].node[1]] += resb[(m-1)*dims+1]; */
+/* break; */
/* case 3: */
-/* RL[j][E->sien[j][els].node[2]] += resb[(m-1)*dims+1]; */
-/* break; */
+/* RL[j][E->sien[j][els].node[2]] += resb[(m-1)*dims+1]; */
+/* break; */
/* case 7: */
-/* RL[j][E->sien[j][els].node[3]] += resb[(m-1)*dims+1]; */
-/* break; */
+/* RL[j][E->sien[j][els].node[3]] += resb[(m-1)*dims+1]; */
+/* break; */
/* case 6: */
-/* RL[j][E->sien[j][els].node[4]] += resb[(m-1)*dims+1]; */
-/* break; */
-/* } */
+/* RL[j][E->sien[j][els].node[4]] += resb[(m-1)*dims+1]; */
+/* break; */
+/* } */
/* } */
@@ -565,7 +780,7 @@
/* dGammax.vpt[GMVGAMMA(1,n)] * l_1d[n].weight[dims-1] */
/* * E->L.vpt[GMVINDEX(m,n)] * E->L.vpt[GMVINDEX(m,n)]; */
/* eltTL[m-1] += */
-/* dGammabx.vpt[GMVGAMMA(1+dims,n)]*l_1d[n].weight[dims-1] */
+/* dGammabx.vpt[GMVGAMMA(1+dims,n)]*l_1d[n].weight[dims-1] */
/* * E->L.vpt[GMVINDEX(m,n)] * E->L.vpt[GMVINDEX(m,n)]; */
/* } */
/* } */
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2006-11-15 22:01:35 UTC (rev 5295)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2006-11-16 01:52:25 UTC (rev 5296)
@@ -232,7 +232,6 @@
struct Shape_function_dA {
double vpt[8+1];
- double spt[4+1];
double ppt[1+1]; };
struct Shape_function1_dA {
@@ -249,7 +248,6 @@
struct Shape_function {
double vpt[8*8]; /* node & gauss pt */
- double spt[8*4]; /* node & gauss pt */
double ppt[8*1]; };
struct Shape_function_dx {
@@ -327,10 +325,14 @@
int total_surf_proc;
int ****loc2proc_map;
+
+#if 0
int nproc_sph[3];
int me_sph;
int me_loc_sph[3];
+#endif
+
int redundant[MAX_LEVELS];
int idb;
int me_loc[4];
@@ -374,6 +376,22 @@
int caps_per_proc;
int capid[NCS];
int max_connections;
+ int *hindex[100];
+ int hindice;
+ float *harm_tpgt[2];
+ float *harm_tpgb[2];
+ float *harm_geoid[2];
+
+ double **tablenplm;
+ double **tablencosf;
+ double **tablensinf;
+
+ double **tablesplm[NCS];
+ double **tablescosf[NCS];
+ double **tablessinf[NCS];
+
+ /**/
+#if 0
int nox;
int noy;
int noz;
@@ -384,22 +402,19 @@
int elz;
int snel;
int llmax;
+
int *int_cap;
int *int_ele;
- int *hindex[100];
- int hindice;
- float *harm_tpgt[2];
- float *harm_tpgb[2];
float *harm_slab[2];
float *harm_velp[2];
float *harm_velt[2];
float *harm_divg[2];
float *harm_vort[2];
float *harm_visc[2];
- float *harm_geoid[2];
double *sx[4];
+ double *det[5];
double *con;
- double *det[5];
+
double *tableplm[181];
double *tablecosf[361];
double *tablesinf[361];
@@ -408,7 +423,9 @@
double *tableplm_n[182];
double *tablecosf_n[362];
double *tablesinf_n[362];
+ /**/
struct SIEN *sien;
+#endif
double area[NCS];
double angle[NCS][5];
@@ -422,7 +439,6 @@
float *radius;
int slab_layers;
- int output_llmax;
int lnox;
@@ -793,6 +809,8 @@
char format[20]; /* ascii-local, ascii or hdf5 */
char optional[1000]; /* comma-delimited list of objects to output */
+ int llmax; /* max degree of spherical harmonics output */
+
/* size of collective buffer used by MPI-IO */
int cb_block_size;
int cb_buffer_size;
@@ -814,6 +832,7 @@
int pressure; /* whether to output pressure */
int surf; /* whether to output surface data */
int botm; /* whether to output bottom data */
+ int geoid; /* whether to output geoid/topo spherial harmonics */
int horiz_avg; /* whether to output horizontal averaged profile */
};
More information about the cig-commits
mailing list