[cig-commits] r8973 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Fri Dec 28 23:02:50 PST 2007
Author: becker
Date: 2007-12-28 23:02:49 -0800 (Fri, 28 Dec 2007)
New Revision: 8973
Modified:
mc/3D/CitcomS/trunk/lib/Construct_arrays.c
mc/3D/CitcomS/trunk/lib/Full_geometry_cartesian.c
mc/3D/CitcomS/trunk/lib/Full_sphere_related.c
mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
mc/3D/CitcomS/trunk/lib/Sphere_util.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Trying to sync back up, those should be old changes, I hope.
Modified: mc/3D/CitcomS/trunk/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -661,7 +661,7 @@
for(j=0;j<E->lmesh.NEQ[lev];j++) {
if(E->BI[lev][m][j] ==0.0) fprintf(stderr,"me= %d level %d, equation %d/%d has zero diagonal term\n",E->parallel.me,lev,j,E->mesh.NEQ[lev]);
assert( E->BI[lev][m][j] != 0 /* diagonal of matrix = 0, not acceptable */);
- E->BI[lev][m][j] = (double) 1.0/E->BI[lev][m][j];
+ E->BI[lev][m][j] = (float) 1.0/E->BI[lev][m][j];
}
} /* end for level */
Modified: mc/3D/CitcomS/trunk/lib/Full_geometry_cartesian.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_geometry_cartesian.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Full_geometry_cartesian.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -85,9 +85,7 @@
E->sphere.caps = 12;
E->sphere.max_connections = 6;
- /* adjust the corner coordinates so that the size (surface area) of
- each cap is about the same. */
- offset = 9.736/180.0*M_PI;
+ offset = 10.0/180.0*M_PI;
for (i=1;i<=4;i++) {
E->sphere.cap[(i-1)*3+1].theta[1] = 0.0;
Modified: mc/3D/CitcomS/trunk/lib/Full_sphere_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_sphere_related.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Full_sphere_related.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -25,8 +25,6 @@
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
-
-
/* Functions relating to the building and use of mesh locations ... */
@@ -35,181 +33,41 @@
#include "element_definitions.h"
#include "global_defs.h"
+void full_coord_of_cap(E,m,icap)
+ struct All_variables *E;
+ int icap,m;
+ {
-/**************************************************************/
-/* This function transforms theta and phi to new coords */
-/* u and v using gnomonic projection. */
-/* See http://mathworld.wolfram.com/GnomonicProjection.html */
-
-void spherical_to_uv2(double center[2], int len,
- double *theta, double *phi,
- double *u, double *v)
-{
- double theta_f, phi_f;
- double cos_tf, sin_tf;
- double cosc, cost, sint, cosp2, sinp2;
- int i;
-
- /* theta_f and phi_f are the reference points at the midpoint of the cap */
-
- theta_f = center[0];
- phi_f = center[1];
-
- cos_tf = cos(theta_f);
- sin_tf = sin(theta_f);
-
- for(i=0; i<len; i++) {
- cost = cos(theta[i]);
- sint = sin(theta[i]);
-
- cosp2 = cos(phi[i] - phi_f);
- sinp2 = sin(phi[i] - phi_f);
-
- cosc = cos_tf * cost + sin_tf * sint * cosp2;
- cosc = 1.0 / cosc;
-
- u[i] = sint * sinp2 * cosc;
- v[i] = (sin_tf * cost - cos_tf * sint * cosp2) * cosc;
- }
- return;
-}
-
-
-/**************************************************************/
-/* This function transforms u and v to spherical coord */
-/* theta and phi using inverse gnomonic projection. */
-/* See http://mathworld.wolfram.com/GnomonicProjection.html */
-
-void uv_to_spherical(double center[2], int len,
- double *u, double *v,
- double *theta, double *phi)
-{
- double theta_f, phi_f, cos_tf, sin_tf;
- double x, y, r, c;
- double cosc, sinc;
- int i;
-
- /* theta_f and phi_f are the reference points at the midpoint of the cap */
-
- theta_f = center[0];
- phi_f = center[1];
-
- cos_tf = cos(theta_f);
- sin_tf = sin(theta_f);
-
- for(i=0; i<len; i++) {
- x = u[i];
- y = v[i];
- r = sqrt(x*x + y*y);
-
- /* special case: r=0, then (u,v) is the reference point */
- if(r == 0) {
- theta[i] = theta_f;
- phi[i] = phi_f;
- continue;
- }
-
- /* c = atan(r); cosc = cos(c); sinc = sin(c);*/
- cosc = 1.0 / sqrt(1 + r*r);
- sinc = sqrt(1 - cosc*cosc);
-
- theta[i] = acos(cosc * cos_tf +
- y * sinc * sin_tf / r);
- phi[i] = phi_f + atan(x * sinc /
- (r * sin_tf * cosc - y * cos_tf * sinc));
- }
- return;
-}
-
-
-/* Find the intersection point of two lines */
-/* The lines are: (x[0], y[0]) to (x[1], y[1]) */
-/* (x[2], y[2]) to (x[3], y[3]) */
-/* If found, the intersection point is stored */
-/* in (px, py) and return 1 */
-/* If not found, return 0 */
-
-static int find_intersection(double *x, double *y,
- double *px, double *py)
-{
- double a1, b1, c1;
- double a2, b2, c2;
- double denom;
-
- a1 = y[1] - y[0];
- b1 = x[0] - x[1];
- c1 = x[1]*y[0] - x[0]*y[1];
-
- a2 = y[3] - y[2];
- b2 = x[2] - x[3];
- c2 = x[3]*y[2] - x[2]*y[3];
-
- denom = a1*b2 - a2*b1;
- if (denom == 0) return 0; /* the lines are parallel! */
-
- *px = (b1*c2 - b2*c1)/denom;
- *py = (a2*c1 - a1*c2)/denom;
- return 1;
-}
-
-
-void full_coord_of_cap(struct All_variables *E, int m, int icap)
-{
- int i, j, k, lev, temp, elx, ely, nox, noy, noz;
- int node, snode;
- int lelx, lely, lnox, lnoy, lnoz;
- int ok;
- double x[5], y[5], z[5], referencep[2];
- double *theta0, *fi0;
- double *tt1, *tt2, *tt3, *tt4, *ff1, *ff2, *ff3, *ff4;
- double *u1, *u2, *u3, *u4, *v1, *v2, *v3, *v4;
- double *px, *py, *qx, *qy;
- double theta, fi, cost, sint, cosf, sinf;
- double a, b;
+ int i,j,k,lev,temp,elx,ely,nox,noy,noz,node,nodes;
+ int lelx,lely,lnox,lnoy,lnoz;
+ double x[5],y[5],z[5],xx[5],yy[5],zz[5];
+ double *theta1,*fi1,*theta2,*fi2,*theta,*fi,*SX[2];
double myatan();
void even_divide_arc12();
temp = max(E->mesh.NOY[E->mesh.levmax],E->mesh.NOX[E->mesh.levmax]);
- theta0 = (double *)malloc((temp+1)*sizeof(double));
- fi0 = (double *)malloc((temp+1)*sizeof(double));
+ theta1 = (double *)malloc((temp+1)*sizeof(double));
+ fi1 = (double *)malloc((temp+1)*sizeof(double));
+ theta2 = (double *)malloc((temp+1)*sizeof(double));
+ fi2 = (double *)malloc((temp+1)*sizeof(double));
+ theta = (double *)malloc((temp+1)*sizeof(double));
+ fi = (double *)malloc((temp+1)*sizeof(double));
- tt1 = (double *)malloc((temp+1)*sizeof(double));
- tt2 = (double *)malloc((temp+1)*sizeof(double));
- tt3 = (double *)malloc((temp+1)*sizeof(double));
- tt4 = (double *)malloc((temp+1)*sizeof(double));
+ temp = E->mesh.NOY[E->mesh.levmax]*E->mesh.NOX[E->mesh.levmax]; /* allocate enough for the entire cap */
- ff1 = (double *)malloc((temp+1)*sizeof(double));
- ff2 = (double *)malloc((temp+1)*sizeof(double));
- ff3 = (double *)malloc((temp+1)*sizeof(double));
- ff4 = (double *)malloc((temp+1)*sizeof(double));
+ SX[0] = (double *)malloc((temp+1)*sizeof(double));
+ SX[1] = (double *)malloc((temp+1)*sizeof(double));
- u1 = (double *)malloc((temp+1)*sizeof(double));
- u2 = (double *)malloc((temp+1)*sizeof(double));
- u3 = (double *)malloc((temp+1)*sizeof(double));
- u4 = (double *)malloc((temp+1)*sizeof(double));
-
- v1 = (double *)malloc((temp+1)*sizeof(double));
- v2 = (double *)malloc((temp+1)*sizeof(double));
- v3 = (double *)malloc((temp+1)*sizeof(double));
- v4 = (double *)malloc((temp+1)*sizeof(double));
-
- temp = E->mesh.NOY[E->mesh.levmax] * E->mesh.NOX[E->mesh.levmax];
- px = malloc((temp+1)*sizeof(double));
- py = malloc((temp+1)*sizeof(double));
- qx = malloc((temp+1)*sizeof(double));
- qy = malloc((temp+1)*sizeof(double));
-
- /* 4 corners of the cap */
for (i=1;i<=4;i++) {
x[i] = sin(E->sphere.cap[icap].theta[i])*cos(E->sphere.cap[icap].fi[i]);
y[i] = sin(E->sphere.cap[icap].theta[i])*sin(E->sphere.cap[icap].fi[i]);
z[i] = cos(E->sphere.cap[icap].theta[i]);
- }
+ }
+
+ for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
-
elx = E->lmesh.ELX[lev]*E->parallel.nprocx;
ely = E->lmesh.ELY[lev]*E->parallel.nprocy;
nox = elx+1;
@@ -220,176 +78,67 @@
lnox = lelx+1;
lnoy = lely+1;
lnoz = E->lmesh.NOZ[lev];
+ /* evenly divide arc linking 1 and 2, and the arc linking 4 and 3 */
- /* evenly divide arc linking 1 and 2 */
- even_divide_arc12(elx,x[1],y[1],z[1],x[2],y[2],z[2],theta0,fi0);
+ even_divide_arc12(elx,x[1],y[1],z[1],x[2],y[2],z[2],theta1,fi1);
+ even_divide_arc12(elx,x[4],y[4],z[4],x[3],y[3],z[3],theta2,fi2);
- /* pick up only points within this processor */
- for (j=0, i=E->lmesh.nxs; j<lnox; j++, i++) {
- tt1[j] = theta0[i];
- ff1[j] = fi0[i];
- }
+ for (j=1;j<=nox;j++) {
- /* evenly divide arc linking 4 and 3 */
- even_divide_arc12(elx,x[4],y[4],z[4],x[3],y[3],z[3],theta0,fi0);
+ /* pick up the two ends */
+ xx[1] = sin(theta1[j])*cos(fi1[j]);
+ yy[1] = sin(theta1[j])*sin(fi1[j]);
+ zz[1] = cos(theta1[j]);
+ xx[2] = sin(theta2[j])*cos(fi2[j]);
+ yy[2] = sin(theta2[j])*sin(fi2[j]);
+ zz[2] = cos(theta2[j]);
- /* pick up only points within this processor */
- for (j=0, i=E->lmesh.nxs; j<lnox; j++, i++) {
- tt2[j] = theta0[i];
- ff2[j] = fi0[i];
- }
- /* evenly divide arc linking 1 and 4 */
- even_divide_arc12(ely,x[1],y[1],z[1],x[4],y[4],z[4],theta0,fi0);
+ even_divide_arc12(ely,xx[1],yy[1],zz[1],xx[2],yy[2],zz[2],theta,fi);
- /* pick up only points within this processor */
- for (k=0, i=E->lmesh.nys; k<lnoy; k++, i++) {
- tt3[k] = theta0[i];
- ff3[k] = fi0[i];
- }
+ for (k=1;k<=noy;k++) {
+ nodes = j + (k-1)*nox;
+ SX[0][nodes] = theta[k];
+ SX[1][nodes] = fi[k];
+ }
- /* evenly divide arc linking 2 and 3 */
- even_divide_arc12(ely,x[2],y[2],z[2],x[3],y[3],z[3],theta0,fi0);
- /* pick up only points within this processor */
- for (k=0, i=E->lmesh.nys; k<lnoy; k++, i++) {
- tt4[k] = theta0[i];
- ff4[k] = fi0[i];
- }
+ } /* end for j */
- /* compute the intersection point of these great circles */
- /* the point is first found in u-v space and project back */
- /* to theta-phi space later */
+ /* get the coordinates for the entire cap */
- referencep[0] = E->sphere.cap[icap].theta[2];
- referencep[1] = E->sphere.cap[icap].fi[2];
+ for (j=1;j<=lnox;j++)
+ for (k=1;k<=lnoy;k++) {
+ nodes = (j+(E->lmesh.NXS[lev]-1))+(k-1+(E->lmesh.NYS[lev]-1))*nox;
+ for (i=1;i<=lnoz;i++) {
+ node = i + (j-1)*lnoz + (k-1)*lnox*lnoz;
- spherical_to_uv2(referencep, lnox, tt1, ff1, u1, v1);
- spherical_to_uv2(referencep, lnox, tt2, ff2, u2, v2);
- spherical_to_uv2(referencep, lnoy, tt3, ff3, u3, v3);
- spherical_to_uv2(referencep, lnoy, tt4, ff4, u4, v4);
+ /* theta,fi,and r coordinates */
+ E->SX[lev][m][1][node] = SX[0][nodes];
+ E->SX[lev][m][2][node] = SX[1][nodes];
+ E->SX[lev][m][3][node] = E->sphere.R[lev][i];
- snode = 0;
- for(k=0; k<lnoy; k++) {
- x[2] = u3[k];
- y[2] = v3[k];
-
- x[3] = u4[k];
- y[3] = v4[k];
-
- for(j=0; j<lnox; j++) {
- x[0] = u1[j];
- y[0] = v1[j];
-
- x[1] = u2[j];
- y[1] = v2[j];
-
- ok = find_intersection(x, y, &a, &b);
- if(!ok) {
- fprintf(stderr, "Error(Full_coord_of_cap): cannot find intersection point!\n");
- exit(10);
+ /* x,y,and z oordinates */
+ E->X[lev][m][1][node] =
+ E->sphere.R[lev][i]*sin(SX[0][nodes])*cos(SX[1][nodes]);
+ E->X[lev][m][2][node] =
+ E->sphere.R[lev][i]*sin(SX[0][nodes])*sin(SX[1][nodes]);
+ E->X[lev][m][3][node] =
+ E->sphere.R[lev][i]*cos(SX[0][nodes]);
+ }
}
- px[snode] = a;
- py[snode] = b;
- snode++;
- }
- }
+ } /* end for lev */
- uv_to_spherical(referencep, snode, px, py, qx, qy);
+ free ((void *)SX[0]);
+ free ((void *)SX[1]);
+ free ((void *)theta );
+ free ((void *)theta1);
+ free ((void *)theta2);
+ free ((void *)fi );
+ free ((void *)fi1 );
+ free ((void *)fi2 );
- /* replace (qx, qy) by (tt?, ff?) for points on the edge */
- if(E->parallel.me_loc[2] == 0) {
- /* x boundary */
- for(k=0; k<lnox; k++) {
- i = k;
- qx[i] = tt1[k];
- qy[i] = ff1[k];
- }
- }
-
- if(E->parallel.me_loc[2] == E->parallel.nprocy-1) {
- /* x boundary */
- for(k=0; k<lnox; k++) {
- i = k + (lnoy - 1) * lnox;
- qx[i] = tt2[k];
- qy[i] = ff2[k];
- }
- }
-
- if(E->parallel.me_loc[1] == 0) {
- /* y boundary */
- for(k=0; k<lnoy; k++) {
- i = k * lnox;
- qx[i] = tt3[k];
- qy[i] = ff3[k];
- }
- }
-
- if(E->parallel.me_loc[1] == E->parallel.nprocx-1) {
- /* y boundary */
- for(k=0; k<lnoy; k++) {
- i = (k + 1) * lnox - 1;
- qx[i] = tt4[k];
- qy[i] = ff4[k];
- }
- }
-
- /* store the node location in spherical and cartesian coordinates */
- for (k=0; k<snode; k++) {
- theta = qx[k];
- fi = qy[k];
-
- cost = cos(theta);
- sint = sin(theta);
- cosf = cos(fi);
- sinf = sin(fi);
-
- for (i=1; i<=lnoz; i++) {
- node = i + k*lnoz;
-
- /* theta,fi,and r coordinates */
- E->SX[lev][m][1][node] = theta;
- E->SX[lev][m][2][node] = fi;
- E->SX[lev][m][3][node] = E->sphere.R[lev][i];
-
- /* x,y,and z oordinates */
- E->X[lev][m][1][node] = E->sphere.R[lev][i]*sint*cosf;
- E->X[lev][m][2][node] = E->sphere.R[lev][i]*sint*sinf;
- E->X[lev][m][3][node] = E->sphere.R[lev][i]*cost;
- }
- }
-
- } /* end for lev */
-
- free ((void *)theta0);
- free ((void *)fi0 );
-
- free ((void *)tt1 );
- free ((void *)tt2 );
- free ((void *)tt3 );
- free ((void *)tt4 );
-
- free ((void *)ff1 );
- free ((void *)ff2 );
- free ((void *)ff3 );
- free ((void *)ff4 );
-
- free ((void *)u1 );
- free ((void *)u2 );
- free ((void *)u3 );
- free ((void *)u4 );
-
- free ((void *)v1 );
- free ((void *)v2 );
- free ((void *)v3 );
- free ((void *)v4 );
-
- free ((void *)px );
- free ((void *)py );
- free ((void *)qx );
- free ((void *)qy );
-
return;
-}
+ }
Modified: mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Full_tracer_advection.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -1834,6 +1834,13 @@
fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
+
+ if ((E->parallel.nprocx > 1) || (E->parallel.nprocy > 1)) {
+ fprintf(E->trace.fpt,"Sorry - Tracer code does not (yet) work if nprocx or nprocy is greater than 1\n");
+ fflush(E->trace.fpt);
+ parallel_process_termination();
+ }
+
if (E->trace.ic_method==0)
{
fprintf(E->trace.fpt,"Generating New Tracer Array\n");
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -126,40 +126,7 @@
}
-
/* =================================================
- rotate the mesh by a rotation matrix
- =================================================*/
-static void full_rotate_mesh(struct All_variables *E, double dircos[4][4],
- int m, int icap)
-{
- int i,lev;
- double t[4], myatan();
-
- for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- for (i=1;i<=E->lmesh.NNO[lev];i++) {
- t[0] = E->X[lev][m][1][i]*dircos[1][1]+
- E->X[lev][m][2][i]*dircos[1][2]+
- E->X[lev][m][3][i]*dircos[1][3];
- t[1] = E->X[lev][m][1][i]*dircos[2][1]+
- E->X[lev][m][2][i]*dircos[2][2]+
- E->X[lev][m][3][i]*dircos[2][3];
- t[2] = E->X[lev][m][1][i]*dircos[3][1]+
- E->X[lev][m][2][i]*dircos[3][2]+
- E->X[lev][m][3][i]*dircos[3][3];
-
- E->X[lev][m][1][i] = t[0];
- E->X[lev][m][2][i] = t[1];
- E->X[lev][m][3][i] = t[2];
- E->SX[lev][m][1][i] = acos(t[2]/E->SX[lev][m][3][i]);
- E->SX[lev][m][2][i] = myatan(t[1],t[0]);
- }
- } /* lev */
-
- return;
-}
-
-/* =================================================
Standard node positions including mesh refinement
================================================= */
@@ -169,7 +136,6 @@
{
int i,j,k,ii,lev;
double ro,dr,*rr,*RR,fo;
- double dircos[4][4];
float tt1;
int step,nn;
char output_file[255], a[255];
@@ -243,54 +209,30 @@
free ((void *) rr);
free ((void *) RR);
+ ro = -0.5*(M_PI/4.0)/E->mesh.elx;
+ fo = 0.0;
+
+ E->sphere.dircos[1][1] = cos(ro)*cos(fo);
+ E->sphere.dircos[1][2] = cos(ro)*sin(fo);
+ E->sphere.dircos[1][3] = -sin(ro);
+ E->sphere.dircos[2][1] = -sin(fo);
+ E->sphere.dircos[2][2] = cos(fo);
+ E->sphere.dircos[2][3] = 0.0;
+ E->sphere.dircos[3][1] = sin(ro)*cos(fo);
+ E->sphere.dircos[3][2] = sin(ro)*sin(fo);
+ E->sphere.dircos[3][3] = cos(ro);
+
for (j=1;j<=E->sphere.caps_per_proc;j++) {
ii = E->sphere.capid[j];
full_coord_of_cap(E,j,ii);
}
-
- if (E->control.verbose) {
- for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- fprintf(E->fp_out,"output_coordinates before rotation %d \n",lev);
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- for (i=1;i<=E->lmesh.NNO[lev];i++)
- if(i%E->lmesh.NOZ[lev]==1)
- fprintf(E->fp_out,"%d %d %g %g %g\n",j,i,E->SX[lev][j][1][i],E->SX[lev][j][2][i],E->SX[lev][j][3][i]);
- }
- fflush(E->fp_out);
- }
-
/* rotate the mesh to avoid two poles on mesh points */
-
- ro = -0.5*(M_PI/4.0)/E->mesh.elx;
- fo = 0.0;
-
- dircos[1][1] = cos(ro)*cos(fo);
- dircos[1][2] = cos(ro)*sin(fo);
- dircos[1][3] = -sin(ro);
- dircos[2][1] = -sin(fo);
- dircos[2][2] = cos(fo);
- dircos[2][3] = 0.0;
- dircos[3][1] = sin(ro)*cos(fo);
- dircos[3][2] = sin(ro)*sin(fo);
- dircos[3][3] = cos(ro);
-
for (j=1;j<=E->sphere.caps_per_proc;j++) {
ii = E->sphere.capid[j];
- full_rotate_mesh(E,dircos,j,ii);
+ rotate_mesh(E,j,ii);
}
- if (E->control.verbose) {
- for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- fprintf(E->fp_out,"output_coordinates after rotation %d \n",lev);
- for (j=1;j<=E->sphere.caps_per_proc;j++)
- for (i=1;i<=E->lmesh.NNO[lev];i++)
- if(i%E->lmesh.NOZ[lev]==1)
- fprintf(E->fp_out,"%d %d %g %g %g\n",j,i,E->SX[lev][j][1][i],E->SX[lev][j][2][i],E->SX[lev][j][3][i]);
- }
- fflush(E->fp_out);
- }
-
compute_angle_surf_area (E); /* used for interpolation */
for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++)
@@ -301,6 +243,22 @@
E->SinCos[lev][j][2][i] = cos(E->SX[lev][j][1][i]);
E->SinCos[lev][j][3][i] = cos(E->SX[lev][j][2][i]);
}
+
+ /*
+if (E->control.verbose) {
+ for (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ fprintf(E->fp_out,"output_coordinates after rotation %d \n",lev);
+ for (j=1;j<=E->sphere.caps_per_proc;j++)
+ for (i=1;i<=E->lmesh.NNO[lev];i++)
+ if(i%E->lmesh.NOZ[lev]==1)
+ fprintf(E->fp_out,"%d %d %g %g %g\n",j,i,E->SX[lev][j][1][i],E->SX[lev][j][2][i],E->SX[lev][j][3][i]);
+ }
+ fflush(E->fp_out);
+}
+ */
+
+
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -1141,7 +1141,10 @@
sprintf(output_file,"%s/qt.dat", E->control.data_dir);
else
sprintf(output_file,"%s.qt.dat", E->control.data_file);
- E->output.fpqt = output_open(output_file, "w");
+ if(E->control.restart)
+ E->output.fpqt = output_open(output_file, "a"); /* append for restart */
+ else
+ E->output.fpqt = output_open(output_file, "w");
}else{
E->output.fpqt = NULL;
}
@@ -1151,7 +1154,10 @@
sprintf(output_file,"%s/qb.dat", E->control.data_dir);
else
sprintf(output_file,"%s.qb.dat", E->control.data_file);
- E->output.fpqb = output_open(output_file, "w");
+ if(E->control.restart)
+ E->output.fpqb = output_open(output_file, "a"); /* append */
+ else
+ E->output.fpqb = output_open(output_file, "w");
}else{
E->output.fpqb = NULL;
}
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -1440,26 +1440,17 @@
status = set_attribute_float(input, "Ra_410", E->control.Ra_410);
status = set_attribute_float(input, "clapeyron410", E->control.clapeyron410);
status = set_attribute_float(input, "transT410", E->control.transT410);
- status = set_attribute_float(input, "width410",
- (E->control.inv_width410 == 0)?
- E->control.inv_width410 :
- 1.0/E->control.inv_width410);
+ status = set_attribute_float(input, "width410", E->control.width410);
status = set_attribute_float(input, "Ra_670", E->control.Ra_670);
status = set_attribute_float(input, "clapeyron670", E->control.clapeyron670);
status = set_attribute_float(input, "transT670", E->control.transT670);
- status = set_attribute_float(input, "width670",
- (E->control.inv_width670 == 0)?
- E->control.inv_width670 :
- 1.0/E->control.inv_width670);
+ status = set_attribute_float(input, "width670", E->control.width670);
status = set_attribute_float(input, "Ra_cmb", E->control.Ra_cmb);
status = set_attribute_float(input, "clapeyroncmb", E->control.clapeyroncmb);
status = set_attribute_float(input, "transTcmb", E->control.transTcmb);
- status = set_attribute_float(input, "widthcmb",
- (E->control.inv_widthcmb == 0)?
- E->control.inv_widthcmb :
- 1.0/E->control.inv_widthcmb);
+ status = set_attribute_float(input, "widthcmb", E->control.widthcmb);
/*
* Solver.inventory
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -229,6 +229,22 @@
} /* lev */
+/* do not need to rotate for the mesh grid for one regional problem */
+
+
+ ro = -0.5*(M_PI/4.0)/E->lmesh.elx;
+ fo = 0.0;
+
+ E->sphere.dircos[1][1] = cos(ro)*cos(fo);
+ E->sphere.dircos[1][2] = cos(ro)*sin(fo);
+ E->sphere.dircos[1][3] = -sin(ro);
+ E->sphere.dircos[2][1] = -sin(fo);
+ E->sphere.dircos[2][2] = cos(fo);
+ E->sphere.dircos[2][3] = 0.0;
+ E->sphere.dircos[3][1] = sin(ro)*cos(fo);
+ E->sphere.dircos[3][2] = sin(ro)*sin(fo);
+ E->sphere.dircos[3][3] = cos(ro);
+
for (j=1;j<=E->sphere.caps_per_proc;j++) {
regional_coord_of_cap(E,j,0);
}
@@ -246,6 +262,12 @@
}
fflush(E->fp_out);
}
+ /* rotate the mesh to avoid two poles on mesh points */
+/*
+ for (j=1;j<=E->sphere.caps_per_proc;j++) {
+ rotate_mesh(E,j,0);
+ }
+*/
compute_angle_surf_area (E); /* used for interpolation */
Modified: mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Sphere_harmonics.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -141,3 +141,195 @@
return;
}
+/* ================================================
+ compute angle and area
+ ================================================*/
+
+void compute_angle_surf_area (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();
+ 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 (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);
+
+ 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);
+
+ 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]];
+ }
+
+ 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];
+
+ 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 */
+
+ } /* end for m */
+
+ return;
+}
+
+/* ================================================
+ area of spherical rectangle
+ ================================================ */
+double area_sphere_cap(angle)
+ double angle[6];
+{
+
+ 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[3];
+ b = angle[4];
+ c = angle[5];
+ area += area_of_sphere_triag(a,b,c);
+
+ return (area);
+}
+
+/* ================================================
+ area of spherical triangle
+ ================================================ */
+double area_of_sphere_triag(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;
+
+ 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);
+}
+
+/* =====================================================================
+ 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];
+
+ for (i=1;i<=4;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];
+ }
+
+ 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]);
+
+ return (area1);
+}
+
+/* ================================
+ get the angle for given four points spherical rectangle
+ ================================= */
+
+void get_angle_sphere_cap(xx,angle)
+ double xx[4][5],angle[6];
+{
+
+ 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 (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;
+}
+
+/* ================================
+ 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);
+}
Modified: mc/3D/CitcomS/trunk/lib/Sphere_util.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Sphere_util.c 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/Sphere_util.c 2007-12-29 07:02:49 UTC (rev 8973)
@@ -60,195 +60,37 @@
return;
}
-/* ================================================
- compute angle and area
- ================================================*/
-void compute_angle_surf_area (E)
- struct All_variables *E;
-{
+/* =================================================
+ rotate the mesh
+ =================================================*/
+void rotate_mesh(E,m,icap)
+ struct All_variables *E;
+ int icap,m;
+ {
- 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();
+ int i,lev;
+ double t[4],myatan();
- 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 (lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ for (i=1;i<=E->lmesh.NNO[lev];i++) {
+ t[0] = E->X[lev][m][1][i]*E->sphere.dircos[1][1]+
+ E->X[lev][m][2][i]*E->sphere.dircos[1][2]+
+ E->X[lev][m][3][i]*E->sphere.dircos[1][3];
+ t[1] = E->X[lev][m][1][i]*E->sphere.dircos[2][1]+
+ E->X[lev][m][2][i]*E->sphere.dircos[2][2]+
+ E->X[lev][m][3][i]*E->sphere.dircos[2][3];
+ t[2] = E->X[lev][m][1][i]*E->sphere.dircos[3][1]+
+ E->X[lev][m][2][i]*E->sphere.dircos[3][2]+
+ E->X[lev][m][3][i]*E->sphere.dircos[3][3];
- 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]];
- }
+ E->X[lev][m][1][i] = t[0];
+ E->X[lev][m][2][i] = t[1];
+ E->X[lev][m][3][i] = t[2];
+ E->SX[lev][m][1][i] = acos(t[2]/E->SX[lev][m][3][i]);
+ E->SX[lev][m][2][i] = myatan(t[1],t[0]);
+ }
+ } /* lev */
- 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];
-
- 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 (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);
-
- 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);
-
-/* 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 */
-
- } /* end for m */
-
- return;
-}
-
-/* ================================================
- area of spherical rectangle
- ================================================ */
-double area_sphere_cap(angle)
- double angle[6];
-{
-
- 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[3];
- b = angle[4];
- c = angle[5];
- area += area_of_sphere_triag(a,b,c);
-
- return (area);
-}
-
-/* ================================================
- area of spherical triangle
- ================================================ */
-double area_of_sphere_triag(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;
-
- 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);
-}
-
-/* =====================================================================
- 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];
-
- for (i=1;i<=4;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];
- }
-
- 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]);
-
- return (area1);
-}
-
-/* ================================
- get the angle for given four points spherical rectangle
- ================================= */
-
-void get_angle_sphere_cap(xx,angle)
- double xx[4][5],angle[6];
-{
-
- 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 (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;
-}
-
-/* ================================
- 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);
-}
+ return;
+ }
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-12-27 02:08:49 UTC (rev 8972)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-12-29 07:02:49 UTC (rev 8973)
@@ -276,6 +276,7 @@
double *area1[MAX_LEVELS][NCS];
double *angle1[MAX_LEVELS][NCS][5];
+ double dircos[4][4];
double *R[MAX_LEVELS];
double ro,ri;
struct CAP cap[NCS];
More information about the cig-commits
mailing list