[cig-commits] r12819 - in mc/3D/CitcomS/branches/v3.0: . lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Sep 5 12:46:12 PDT 2008
Author: tan2
Date: 2008-09-05 12:46:12 -0700 (Fri, 05 Sep 2008)
New Revision: 12819
Modified:
mc/3D/CitcomS/branches/v3.0/
mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c
mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c
mc/3D/CitcomS/branches/v3.0/lib/Size_does_matter.c
mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c
mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c
Log:
Merged revisions 12272-12276 via svnmerge from
svn+ssh://svn@geodynamics.org/cig/mc/3D/CitcomS/trunk
........
r12272 | tan2 | 2008-06-18 18:04:06 -0700 (Wed, 18 Jun 2008) | 6 lines
Compute E->surf_det (top surface area, ie. surface jacobian, of an element) on a unit sphere.
Use E->surf_det, instead of the surface jacobian returned by get_global_1d_shape_fn(), to expand spherical harmonics.
The surface jacobian returned by get_global_1d_shape_fn() is defined at the top radius of the current mesh. When nprocz>1, the top radius is not the same as the outer radius of the sphere for processors. As a result, the spherical harmonics expansion gave incorrect result when nprocz>1 in the previous revisions.
........
r12273 | tan2 | 2008-06-18 18:05:12 -0700 (Wed, 18 Jun 2008) | 2 lines
Send geoid arrays to vertical columns of processors. After this fix, the geoid output is the same for nprocz=1 and nprocz>1.
........
r12276 | tan2 | 2008-06-19 14:42:27 -0700 (Thu, 19 Jun 2008) | 2 lines
Added debugging (but disabled) function for spherical harmonics expansion.
........
Property changes on: mc/3D/CitcomS/branches/v3.0
___________________________________________________________________
Name: svnmerge-integrated
- /mc/3D/CitcomS/trunk:1-8274,8276-8280,8893,9216,11219,11280,11332,12062-12176,12178-12201,12203,12207,12210-12256,12258,12260-12269,12271,12274-12275,12278-12279
+ /mc/3D/CitcomS/trunk:1-8274,8276-8280,8893,9216,11219,11280,11332,12062-12176,12178-12201,12203,12207,12210-12256,12258,12260-12269,12271-12276,12278-12279
Modified: mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c 2008-09-05 17:20:07 UTC (rev 12818)
+++ mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c 2008-09-05 19:46:12 UTC (rev 12819)
@@ -273,6 +273,7 @@
double con;
double modified_plgndr_a(int, int, double);
void temperatures_conform_bcs();
+ void debug_sphere_expansion(struct All_variables *);
noy=E->lmesh.noy;
nox=E->lmesh.nox;
@@ -388,6 +389,8 @@
temperatures_conform_bcs(E);
+ /* debugging the code of expanding spherical harmonics */
+ debug_sphere_expansion(E);
return;
}
Modified: mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c 2008-09-05 17:20:07 UTC (rev 12818)
+++ mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c 2008-09-05 19:46:12 UTC (rev 12819)
@@ -787,6 +787,46 @@
return;
}
+
+/* ================================================== */
+/* ================================================== */
+void broadcast_vertical(struct All_variables *E,
+ float *sphc, float *sphs,
+ int root)
+{
+ int jumpp, total, j;
+ float *temp;
+
+ if(E->parallel.nprocz == 1) return;
+
+ jumpp = E->sphere.hindice;
+ total = E->sphere.hindice*2;
+ temp = (float *) malloc(total*sizeof(float));
+
+ if (E->parallel.me_loc[3] == root) {
+ /* pack */
+ for (j=0; j<E->sphere.hindice; j++) {
+ temp[j] = sphc[j];
+ temp[j+jumpp] = sphs[j];
+ }
+ }
+
+ MPI_Bcast(temp, total, MPI_FLOAT, root, E->parallel.vertical_comm);
+
+ if (E->parallel.me_loc[3] != root) {
+ /* unpack */
+ for (j=0; j<E->sphere.hindice; j++) {
+ sphc[j] = temp[j];
+ sphs[j] = temp[j+jumpp];
+ }
+ }
+
+ free((void *)temp);
+
+ return;
+}
+
+
/*
* remove rigid body rotation from the velocity
*/
Modified: mc/3D/CitcomS/branches/v3.0/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Size_does_matter.c 2008-09-05 17:20:07 UTC (rev 12818)
+++ mc/3D/CitcomS/branches/v3.0/lib/Size_does_matter.c 2008-09-05 19:46:12 UTC (rev 12819)
@@ -275,14 +275,16 @@
const int oned = onedvpoints[E->mesh.nsd];
- double xx[4][5],dxda[4][4];
+ double xx[4][5], dxda[4][4], r2;
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(k=1;k<=oned;k++) { /* all of the vpoints*/
E->surf_det[m][k] = (double *)malloc((1+E->lmesh.snel)*sizeof(double));
}
- for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ r2 = 1.0 / (E->sx[m][3][E->lmesh.elz+1] * E->sx[m][3][E->lmesh.elz+1]);
+
for (es=1;es<=E->lmesh.snel;es++) {
el = es * E->lmesh.elz;
get_side_x_cart(E, xx, el, SIDE_TOP, m);
@@ -298,10 +300,12 @@
dxda[d][e] += xx[e][i]*E->Mx.vpt[GMVXINDEX(d-1,i,k)];
jacobian = determinant(dxda,E->mesh.nsd-1);
- E->surf_det[m][k][es] = jacobian;
+
+ /* scale the jacobian so that it is defined on a unit sphere */
+ E->surf_det[m][k][es] = jacobian * r2;
}
}
-
+ }
return;
}
Modified: mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c 2008-09-05 17:20:07 UTC (rev 12818)
+++ mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c 2008-09-05 19:46:12 UTC (rev 12819)
@@ -49,11 +49,7 @@
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;
for (i=0;i<E->sphere.hindice;i++) {
sphc[i] = 0.0;
@@ -62,10 +58,7 @@
for (m=1;m<=E->sphere.caps_per_proc;m++)
for (es=1;es<=E->lmesh.snel;es++) {
- el = es*E->lmesh.elz;
- 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++) {
@@ -78,12 +71,12 @@
* E->sphere.tablesplm[m][j][p]
* E->sphere.tablescosf[m][j][mm]
* E->M.vpt[GMVINDEX(d,nint)]
- * dGamma.vpt[GMVGAMMA(1,nint)];
+ * E->surf_det[m][nint][es];
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)];
+ * E->surf_det[m][nint][es];
}
}
@@ -97,6 +90,51 @@
}
+void debug_sphere_expansion(struct All_variables *E)
+{
+ /* expand temperature field (which should be a sph. harm. load)
+ * and output the expansion coeff. to stderr
+ */
+ int m, i, j, k, p, node;
+ int ll, mm;
+ float *TT[NCS], *sph_harm[2];
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ TT[m] = (float *) malloc ((E->lmesh.nsf+1)*sizeof(float));
+
+ /* sin coeff */
+ sph_harm[0] = (float*)malloc(E->sphere.hindice*sizeof(float));
+ /* cos coeff */
+ sph_harm[1] = (float*)malloc(E->sphere.hindice*sizeof(float));
+
+ 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->T[m][node];
+ }
+
+ /* expand TT into spherical harmonics */
+ sphere_expansion(E, TT, sph_harm[0], sph_harm[1]);
+
+ /* only the first nprocz CPU needs output */
+ if(E->parallel.me < E->parallel.nprocz) {
+ for (ll=0;ll<=E->output.llmax;ll++)
+ for (mm=0; mm<=ll; mm++) {
+ p = E->sphere.hindex[ll][mm];
+ fprintf(stderr, "T expanded layer=%d ll=%d mm=%d -- %12g %12g\n",
+ k+E->lmesh.nzs-1, ll, mm,
+ sph_harm[0][p], sph_harm[1][p]);
+ }
+ }
+ }
+
+ return;
+}
+
+
/* ==================================================*/
/* ==================================================*/
static void compute_sphereh_table(E)
Modified: mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c 2008-09-05 17:20:07 UTC (rev 12818)
+++ mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c 2008-09-05 19:46:12 UTC (rev 12819)
@@ -30,6 +30,12 @@
#include "element_definitions.h"
#include "global_defs.h"
+void myerror(char *,struct All_variables *);
+void sphere_expansion(struct All_variables *, float **, float *, float *);
+void sum_across_depth_sph1(struct All_variables *, float *, float *);
+void broadcast_vertical(struct All_variables *, float *, float *, int);
+
+
void get_STD_topo(E,tpg,tpgb,divg,vort,ii)
struct All_variables *E;
float **tpg,**tpgb;
@@ -616,11 +622,10 @@
}
}
+ /* send arrays to all processors in the same vertical column */
+ broadcast_vertical(E, geoid_tpgb[0], geoid_tpgb[1], 0);
+ broadcast_vertical(E, geoid_tpgt[0], geoid_tpgt[1], E->parallel.nprocz-1);
- /* accumulate geoid to the surface (top processors) */
- sum_across_depth_sph1(E, geoid_tpgb[0], geoid_tpgb[1]);
-
-
for (j=0; j<2; j++) {
free ((void *) tpgt[j]);
free ((void *) tpgb[j]);
More information about the cig-commits
mailing list