[cig-commits] r9078 - mc/3D/CitcomS/trunk/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Jan 16 12:30:09 PST 2008
Author: tan2
Date: 2008-01-16 12:30:09 -0800 (Wed, 16 Jan 2008)
New Revision: 9078
Modified:
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/lib/tracer_defs.h
Log:
A more memory efficient way to find shape function for tracers.
- Migrated E->tracer.UV to E->gnomonic. The reference point is the 1st node of
the local mesh. The 0th element of E->gnomonic stores the sine of cosine of
reference theta.
- Convertd shape_coefs array from size of nel to size of snel to save memory.
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2008-01-16 19:32:32 UTC (rev 9077)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2008-01-16 20:30:09 UTC (rev 9078)
@@ -837,8 +837,7 @@
}
-/******** GET VELOCITY ***************************************/
-/********************** GNOMONIC INTERPOLATION *******************************/
+/************************ GET VELOCITY ***************************************/
/* */
/* This function interpolates tracer velocity using gnominic interpolation. */
/* Real theta,phi,rad space is transformed into u,v space. This transformation */
@@ -851,7 +850,7 @@
/* to obtain). It was tested that nodal configuration is indeed transformed */
/* into straight lines. */
/* The transformations require a reference point along each cap. Here, the */
-/* midpoint is used. */
+/* first point is used. */
/* Radial and azimuthal shape functions are decoupled. First find the shape */
/* functions associated with the 2D surface plane, then apply radial shape */
/* functions. */
@@ -878,7 +877,7 @@
double *velocity_vector)
{
int iwedge,inum;
- int kk;
+ int i, kk;
int ival;
int itry;
int sphere_key = 0;
@@ -938,8 +937,10 @@
fprintf(E->trace.fpt,"shape %f %f %f\n",shape2d[1],shape2d[2],shape2d[3]);
fprintf(E->trace.fpt,"u %f v %f element: %d \n",u,v, nelem);
fprintf(E->trace.fpt,"Element uv boundaries: \n");
- for(kk=1;kk<=4;kk++)
- fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,E->trace.UV[j][1][E->ien[j][nelem].node[kk]],E->trace.UV[j][2][E->ien[j][nelem].node[kk]]);
+ for(kk=1;kk<=4;kk++) {
+ i = (E->ien[j][nelem].node[kk] - 1) / E->lmesh.noz + 1;
+ fprintf(E->trace.fpt,"%d: U: %f V:%f\n",kk,E->gnomonic[i].u,E->gnomonic[i].v);
+ }
fprintf(E->trace.fpt,"theta: %f phi: %f rad: %f\n",theta,phi,rad);
fprintf(E->trace.fpt,"Element theta-phi boundaries: \n");
for(kk=1;kk<=4;kk++)
@@ -1058,28 +1059,30 @@
{
double a0,a1,a2;
+ /* convert nelem to surface element number */
+ int n = (nelem - 1) / E->lmesh.elz + 1;
/* shape function 1 */
- a0=E->trace.shape_coefs[j][iwedge][1][nelem];
- a1=E->trace.shape_coefs[j][iwedge][2][nelem];
- a2=E->trace.shape_coefs[j][iwedge][3][nelem];
+ a0=E->trace.shape_coefs[j][iwedge][1][n];
+ a1=E->trace.shape_coefs[j][iwedge][2][n];
+ a2=E->trace.shape_coefs[j][iwedge][3][n];
shape2d[1]=a0+a1*u+a2*v;
/* shape function 2 */
- a0=E->trace.shape_coefs[j][iwedge][4][nelem];
- a1=E->trace.shape_coefs[j][iwedge][5][nelem];
- a2=E->trace.shape_coefs[j][iwedge][6][nelem];
+ a0=E->trace.shape_coefs[j][iwedge][4][n];
+ a1=E->trace.shape_coefs[j][iwedge][5][n];
+ a2=E->trace.shape_coefs[j][iwedge][6][n];
shape2d[2]=a0+a1*u+a2*v;
/* shape function 3 */
- a0=E->trace.shape_coefs[j][iwedge][7][nelem];
- a1=E->trace.shape_coefs[j][iwedge][8][nelem];
- a2=E->trace.shape_coefs[j][iwedge][9][nelem];
+ a0=E->trace.shape_coefs[j][iwedge][7][n];
+ a1=E->trace.shape_coefs[j][iwedge][8][n];
+ a2=E->trace.shape_coefs[j][iwedge][9][n];
shape2d[3]=a0+a1*u+a2*v;
@@ -1103,8 +1106,8 @@
double top_bound=1.0+eps;
double bottom_bound=0.0-eps;
- node1=E->IEN[E->mesh.levmax][j][nelem].node[1];
- node5=E->IEN[E->mesh.levmax][j][nelem].node[5];
+ node1=E->ien[j][nelem].node[1];
+ node5=E->ien[j][nelem].node[5];
rad1=E->sx[j][3][node1];
rad5=E->sx[j][3][node5];
@@ -1160,13 +1163,12 @@
double cos_theta_f,sin_theta_f;
double cost,sint,cosp2,sinp2;
- /* theta_f and phi_f are the reference points at the midpoint of the cap */
+ /* theta_f and phi_f are the reference points at the 1st node of the cap */
- theta_f=E->trace.UV[j][1][0];
- phi_f=E->trace.UV[j][2][0];
+ phi_f = E->sx[j][2][1];
- cos_theta_f=E->trace.cos_theta_f;
- sin_theta_f=E->trace.sin_theta_f;
+ cos_theta_f = E->gnomonic[0].u;
+ sin_theta_f = E->gnomonic[0].v;
cost=cos(theta);
/*
@@ -2775,78 +2777,57 @@
/* This function defines nodal points as orthodrome coordinates */
/* u and v. In uv space, great circles form straight lines. */
/* This is used for interpolation method 1 */
-/* UV[j][1][node]=u */
-/* UV[j][2][node]=v */
+/* E->gnomonic[node].u = u */
+/* E->gnomonic[node].v = v */
static void define_uv_space(struct All_variables *E)
{
+ const int j = 1;
+ const int lev = E->mesh.levmax;
+ const int refnode = 1;
+ int i, n;
- int kk,j;
- int midnode;
- int numnodes,node;
+ double u, v, cosc, theta_f, phi_f, dphi, cosd;
+ double *cost, *sint, *cosf, *sinf;
- double u,v,cosc,theta,phi;
- double theta_f,phi_f;
+ if ((E->gnomonic = malloc((E->lmesh.nsf+1)*sizeof(struct GNOMONIC)))
+ == NULL) {
+ fprintf(stderr,"Error(define uv)-not enough memory(a)\n");
+ exit(10);
+ }
- if (E->parallel.me==0) fprintf(stderr,"Setting up UV space\n");
+ sint = E->SinCos[lev][j][0];
+ sinf = E->SinCos[lev][j][1];
+ cost = E->SinCos[lev][j][2];
+ cosf = E->SinCos[lev][j][3];
- numnodes=E->lmesh.nno;
+ /* uv space requires a reference point, using the 1st node */
- /* open memory for uv space coords */
+ theta_f = E->sx[j][1][refnode];
+ phi_f = E->sx[j][2][refnode];
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ /* store cos(theta_f) and sin(theta_f) */
+ E->gnomonic[0].u = cost[refnode];
+ E->gnomonic[0].v = sint[refnode];
- for (kk=1;kk<=2;kk++)
- {
- //TODO: allocate for surface nodes only to save memory
- if ((E->trace.UV[j][kk]=(double *)malloc((numnodes+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"Error(define uv)-not enough memory(a)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
+ /* convert each nodal point to u and v */
- /* uv space requires a reference point */
- /* UV[j][1][0]=fixed theta */
- /* UV[j][2][0]=fixed phi */
+ for (i=1, n=1; i<=E->lmesh.nsf; i++, n+=E->lmesh.noz) {
+ dphi = E->sx[j][2][n] - phi_f;
+ cosd = cos(dphi);
+ cosc = cost[refnode] * cost[n] + sint[refnode] * sint[n] * cosd;
+ u = sint[n] * sin(dphi) / cosc;
+ v = (sint[refnode] * cost[n] - cost[refnode] * sint[n] * cosd)
+ / cosc;
+ E->gnomonic[i].u = u;
+ E->gnomonic[i].v = v;
+ }
- midnode=numnodes/2;
-
- E->trace.UV[j][1][0]=E->sx[j][1][midnode];
- E->trace.UV[j][2][0]=E->sx[j][2][midnode];
-
- theta_f=E->sx[j][1][midnode];
- phi_f=E->sx[j][2][midnode];
-
- E->trace.cos_theta_f=cos(theta_f);
- E->trace.sin_theta_f=sin(theta_f);
-
- /* convert each nodal point to u and v */
-
- for (node=1;node<=numnodes;node++)
- {
- theta=E->sx[j][1][node];
- phi=E->sx[j][2][node];
-
- cosc=cos(theta_f)*cos(theta)+sin(theta_f)*sin(theta)*
- cos(phi-phi_f);
- u=sin(theta)*sin(phi-phi_f)/cosc;
- v=(sin(theta_f)*cos(theta)-cos(theta_f)*sin(theta)*cos(phi-phi_f))/cosc;
-
- E->trace.UV[j][1][node]=u;
- E->trace.UV[j][2][node]=v;
-
- }
-
-
- }/*end cap */
-
return;
}
+
/***************************************************************/
/* DETERMINE SHAPE COEFFICIENTS */
/* */
@@ -2860,118 +2841,97 @@
static void determine_shape_coefficients(struct All_variables *E)
{
+ const int j = 1;
+ int nelem, iwedge, kk, i;
+ int snode;
- int j,nelem,iwedge,kk;
- int node;
+ double u[5], v[5];
+ double x1 = 0.0;
+ double x2 = 0.0;
+ double x3 = 0.0;
+ double y1 = 0.0;
+ double y2 = 0.0;
+ double y3 = 0.0;
+ double delta, a0, a1, a2;
- double u[5],v[5];
- double x1=0.0;
- double x2=0.0;
- double x3=0.0;
- double y1=0.0;
- double y2=0.0;
- double y3=0.0;
- double delta,a0,a1,a2;
+ /* first, allocate memory */
- /* really only have to do this for each surface element, but */
- /* for simplicity, it is done for every element */
+ for(iwedge=1; iwedge<=2; iwedge++) {
+ for (kk=1; kk<=9; kk++) {
+ if ((E->trace.shape_coefs[j][iwedge][kk] =
+ (double *)malloc((E->lmesh.snel+1)*sizeof(double))) == NULL) {
+ fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
+ fflush(E->trace.fpt);
+ exit(10);
+ }
+ }
+ }
- if (E->parallel.me==0) fprintf(stderr," Determining Shape Coefficients\n");
+ for (i=1, nelem=1; i<=E->lmesh.snel; i++, nelem+=E->lmesh.elz) {
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- {
+ /* find u,v of local nodes at one radius */
- /* first, allocate memory */
+ for(kk=1; kk<=4; kk++) {
+ snode = (E->ien[j][nelem].node[kk]-1) / E->lmesh.noz + 1;
+ u[kk] = E->gnomonic[snode].u;
+ v[kk] = E->gnomonic[snode].v;
+ }
- for(iwedge=1;iwedge<=2;iwedge++)
- {
- for (kk=1;kk<=9;kk++)
- {
- //TODO: allocate for surface elements only to save memory
- if ((E->trace.shape_coefs[j][iwedge][kk]=
- (double *)malloc((E->lmesh.nel+1)*sizeof(double)))==NULL)
- {
- fprintf(E->trace.fpt,"ERROR(find shape coefs)-not enough memory(a)\n");
- fflush(E->trace.fpt);
- exit(10);
- }
- }
- }
+ for(iwedge=1; iwedge<=2; iwedge++) {
- for (nelem=1;nelem<=E->lmesh.nel;nelem++)
- {
+ if (iwedge == 1) {
+ x1 = u[1];
+ x2 = u[2];
+ x3 = u[3];
+ y1 = v[1];
+ y2 = v[2];
+ y3 = v[3];
+ }
+ if (iwedge == 2) {
+ x1 = u[1];
+ x2 = u[3];
+ x3 = u[4];
+ y1 = v[1];
+ y2 = v[3];
+ y3 = v[4];
+ }
- /* find u,v of local nodes at one radius */
+ /* shape function 1 */
- for(kk=1;kk<=4;kk++)
- {
- node=E->IEN[E->mesh.levmax][j][nelem].node[kk];
- u[kk]=E->trace.UV[j][1][node];
- v[kk]=E->trace.UV[j][2][node];
- }
+ delta = (x3-x2)*(y1-y2)-(y2-y3)*(x2-x1);
+ a0 = (x2*y3-x3*y2)/delta;
+ a1 = (y2-y3)/delta;
+ a2 = (x3-x2)/delta;
- for(iwedge=1;iwedge<=2;iwedge++)
- {
+ E->trace.shape_coefs[j][iwedge][1][i] = a0;
+ E->trace.shape_coefs[j][iwedge][2][i] = a1;
+ E->trace.shape_coefs[j][iwedge][3][i] = a2;
+ /* shape function 2 */
- if (iwedge==1)
- {
- x1=u[1];
- x2=u[2];
- x3=u[3];
- y1=v[1];
- y2=v[2];
- y3=v[3];
- }
- if (iwedge==2)
- {
- x1=u[1];
- x2=u[3];
- x3=u[4];
- y1=v[1];
- y2=v[3];
- y3=v[4];
- }
+ delta = (x3-x1)*(y2-y1)-(y1-y3)*(x1-x2);
+ a0 = (x1*y3-x3*y1)/delta;
+ a1 = (y1-y3)/delta;
+ a2 = (x3-x1)/delta;
- /* shape function 1 */
+ E->trace.shape_coefs[j][iwedge][4][i] = a0;
+ E->trace.shape_coefs[j][iwedge][5][i] = a1;
+ E->trace.shape_coefs[j][iwedge][6][i] = a2;
- delta=(x3-x2)*(y1-y2)-(y2-y3)*(x2-x1);
- a0=(x2*y3-x3*y2)/delta;
- a1=(y2-y3)/delta;
- a2=(x3-x2)/delta;
+ /* shape function 3 */
- E->trace.shape_coefs[j][iwedge][1][nelem]=a0;
- E->trace.shape_coefs[j][iwedge][2][nelem]=a1;
- E->trace.shape_coefs[j][iwedge][3][nelem]=a2;
+ delta = (x1-x2)*(y3-y2)-(y2-y1)*(x2-x3);
+ a0 = (x2*y1-x1*y2)/delta;
+ a1 = (y2-y1)/delta;
+ a2 = (x1-x2)/delta;
- /* shape function 2 */
+ E->trace.shape_coefs[j][iwedge][7][i] = a0;
+ E->trace.shape_coefs[j][iwedge][8][i] = a1;
+ E->trace.shape_coefs[j][iwedge][9][i] = a2;
- delta=(x3-x1)*(y2-y1)-(y1-y3)*(x1-x2);
- a0=(x1*y3-x3*y1)/delta;
- a1=(y1-y3)/delta;
- a2=(x3-x1)/delta;
+ } /* end wedge */
+ } /* end elem */
- E->trace.shape_coefs[j][iwedge][4][nelem]=a0;
- E->trace.shape_coefs[j][iwedge][5][nelem]=a1;
- E->trace.shape_coefs[j][iwedge][6][nelem]=a2;
-
- /* shape function 3 */
-
- delta=(x1-x2)*(y3-y2)-(y2-y1)*(x2-x3);
- a0=(x2*y1-x1*y2)/delta;
- a1=(y2-y1)/delta;
- a2=(x1-x2)/delta;
-
- E->trace.shape_coefs[j][iwedge][7][nelem]=a0;
- E->trace.shape_coefs[j][iwedge][8][nelem]=a1;
- E->trace.shape_coefs[j][iwedge][9][nelem]=a2;
-
-
- } /* end wedge */
- } /* end elem */
- } /* end cap */
-
-
return;
}
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2008-01-16 19:32:32 UTC (rev 9077)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2008-01-16 20:30:09 UTC (rev 9078)
@@ -654,6 +654,13 @@
};
+struct GNOMONIC {
+ /* gnomonic projected coordinate */
+ double u;
+ double v;
+};
+
+
#ifdef USE_HDF5
#include "hdf5_related.h"
#endif
@@ -691,6 +698,8 @@
/* for chemical convection & composition rheology */
struct COMPOSITION composition;
+ struct GNOMONIC *gnomonic;
+
struct COORD *eco[NCS];
struct IEN *ien[NCS]; /* global */
struct SIEN *sien[NCS];
Modified: mc/3D/CitcomS/trunk/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/tracer_defs.h 2008-01-16 19:32:32 UTC (rev 9077)
+++ mc/3D/CitcomS/trunk/lib/tracer_defs.h 2008-01-16 20:30:09 UTC (rev 9078)
@@ -116,9 +116,6 @@
int *regtoel[13][5];
/* gnomonic shape functions */
- double *UV[13][3];
- double cos_theta_f;
- double sin_theta_f;
double *shape_coefs[13][3][10];
More information about the cig-commits
mailing list