[cig-commits] r7864 - mc/3D/CitcomS/trunk/lib
becker at geodynamics.org
becker at geodynamics.org
Tue Aug 21 18:00:12 PDT 2007
Author: becker
Date: 2007-08-21 18:00:12 -0700 (Tue, 21 Aug 2007)
New Revision: 7864
Added:
mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c
Modified:
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Makefile.am
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_gzdir.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Tracer_setup.c
mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Added an option to remove net rotations of the whole model before
ascii-gz output. This might be something of interest for other output
options as well, but is still experimental. Obvisouly, NR motions
could be removed in a post-processing step, but it is useful to
monitor how those develop during a model run. (Particularly pronounced
for strong lateral viscosity variations based on composition.)
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-22 01:00:12 UTC (rev 7864)
@@ -517,6 +517,7 @@
/* should we add a compositional contribution? */
if(E->control.tracer_enriched){
/* XXX: change Q and Q0 to be a vector of ncomp elements */
+
/* Q = Q0 for C = 0, Q = Q0ER for C = 1, and linearly in
between */
Q *= (1.0 - E->composition.comp_el[m][0][el]);
Added: mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Determine_net_rotation.c 2007-08-22 01:00:12 UTC (rev 7864)
@@ -0,0 +1,447 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<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>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+/*
+
+routines to determine the net rotation velocity of the whole model
+
+TWB
+
+
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "element_definitions.h"
+#include "global_defs.h"
+#include "parsing.h"
+#include "output.h"
+
+double determine_netr_tp(float, float, float, float, float, int, double *, double *);
+void sub_netr(float, float, float, float *, float *, double *);
+void hc_ludcmp_3x3(double [3][3], int *);
+void hc_lubksb_3x3(double [3][3], int *, double *);
+void xyz2rtp(float ,float ,float ,float *);
+void *safe_malloc (size_t );
+double determine_model_net_rotation(struct All_variables *,double *);
+
+/*
+
+determine the mean net rotation of the velocities at all layers
+
+
+modeled after horizontal layer average routines
+
+*/
+double determine_model_net_rotation(struct All_variables *E,double *omega)
+{
+ const int dims = E->mesh.nsd;
+ int m,i,j,k,d,nint,el,elz,elx,ely,j1,j2,i1,i2,k1,k2,nproc;
+ int top,lnode[5],elz9;
+ double *acoef,*coef,*lomega,ddummy,oamp,lamp;
+ float v[2],vw,vtmp,r,t,p,x[3],xp[3],r1,r2,rr;
+ struct Shape_function1 M;
+ struct Shape_function1_dA dGamma;
+ void get_global_1d_shape_fn();
+
+ elz = E->lmesh.elz;elx = E->lmesh.elx;ely = E->lmesh.ely;
+
+ elz9 = elz*9;
+
+ acoef = (double *)safe_malloc(elz9*sizeof(double));
+ coef = (double *)safe_malloc(elz9*sizeof(double));
+ lomega = (double *)safe_malloc(3*elz*sizeof(double));
+
+ for (i=1;i <= elz;i++) { /* loop through depths */
+
+ /* zero out coef for init */
+ determine_netr_tp(ddummy,ddummy,ddummy,ddummy,ddummy,0,(coef+(i-1)*9),&ddummy);
+
+ if (i==elz)
+ top = 1;
+ else
+ top = 0;
+ for (m=1;m <= E->sphere.caps_per_proc;m++)
+ for (k=1;k <= ely;k++)
+ for (j=1;j <= elx;j++) {
+ el = i + (j-1)*elz + (k-1)*elx*elz;
+ get_global_1d_shape_fn(E,el,&M,&dGamma,top,m);
+
+ /* find mean element location and horizontal velocity */
+
+ x[0] = x[1] = x[2] = v[0] = v[1] = vw = 0.0;
+
+ lnode[1] = E->ien[m][el].node[1];
+ lnode[2] = E->ien[m][el].node[2];
+ lnode[3] = E->ien[m][el].node[3];
+ lnode[4] = E->ien[m][el].node[4];
+
+ for(nint=1;nint <= onedvpoints[E->mesh.nsd];nint++) {
+ for(d=1;d <= onedvpoints[E->mesh.nsd];d++){
+ vtmp = E->M.vpt[GMVINDEX(d,nint)] * dGamma.vpt[GMVGAMMA(0,nint)];
+ x[0] += E->x[m][1][lnode[d]] * vtmp; /* coords */
+ x[1] += E->x[m][2][lnode[d]] * vtmp;
+ x[2] += E->x[m][3][lnode[d]] * vtmp;
+ v[0] += E->sphere.cap[m].V[1][lnode[d]] * vtmp; /* theta */
+ v[1] += E->sphere.cap[m].V[2][lnode[d]] * vtmp; /* phi */
+ vw += dGamma.vpt[GMVGAMMA(0,nint)];
+ }
+ }
+ if (i==elz) {
+ lnode[1] = E->ien[m][el].node[5];
+ lnode[2] = E->ien[m][el].node[6];
+ lnode[3] = E->ien[m][el].node[7];
+ lnode[4] = E->ien[m][el].node[8];
+
+ for(nint=1;nint<=onedvpoints[E->mesh.nsd];nint++) {
+ for(d=1;d<=onedvpoints[E->mesh.nsd];d++){
+ vtmp = E->M.vpt[GMVINDEX(d,nint)] * dGamma.vpt[GMVGAMMA(1,nint)];
+ x[0] += E->x[m][1][lnode[d]] * vtmp; /* coords */
+ x[1] += E->x[m][2][lnode[d]] * vtmp;
+ x[2] += E->x[m][3][lnode[d]] * vtmp;
+ v[0] += E->sphere.cap[m].V[1][lnode[d]] * vtmp;
+ v[1] += E->sphere.cap[m].V[2][lnode[d]] * vtmp;
+ vw += dGamma.vpt[GMVGAMMA(1,nint)];
+ }
+ }
+ } /* end of if i==elz */
+ x[0] /= vw;x[1] /= vw;x[2] /= vw; /* convert */
+ xyz2rtp(x[0],x[1],x[2],xp);
+ v[0] /= vw;v[1] /= vw;
+ /* add */
+ determine_netr_tp(xp[0],xp[1],xp[2],v[0],v[1],1,(coef+(i-1)*9),&ddummy);
+ } /* end of j and k, and m */
+
+ } /* Done for i */
+ /*
+ sum it all up
+ */
+ MPI_Allreduce(coef,acoef,elz9,MPI_DOUBLE,MPI_SUM,E->parallel.horizontal_comm);
+
+ omega[0]=omega[1]=omega[2]=0.0;
+
+ /* depth range */
+ rr = E->sx[1][3][E->ien[1][elz].node[5]] - E->sx[1][3][E->ien[1][1].node[1]];
+ vw = 0.0;
+ for (i=0;i < elz;i++) { /* regular 0..n-1 loop */
+ /* solve layer NR */
+ lamp = determine_netr_tp(ddummy,ddummy,ddummy,ddummy,ddummy,2,(acoef+i*9),(lomega+i*3));
+ r1 = E->sx[1][3][E->ien[1][i+1].node[1]]; /* nodal radii for the
+ i-th element, this
+ assumes that there
+ are no lateral
+ variations in radii!
+ */
+ r2 = E->sx[1][3][E->ien[1][i+1].node[5]];
+ vtmp = (r2-r1)/rr; /* weight for this layer */
+ //if(E->parallel.me == 0)
+ // fprintf(stderr,"NR layer %5i (%11g - %11g, %11g): |%11g %11g %11g| = %11g\n",
+ // i+1,r1,r2,vtmp,lomega[i*3+0],lomega[i*3+1],lomega[i*3+2],lamp);
+ /* */
+ for(i1=0;i1<3;i1++)
+ omega[i1] += lomega[i*3+i1] * vtmp;
+ vw += vtmp;
+ }
+ for(i1=0;i1 < 3;i1++)
+ omega[i1] /= vw;
+
+ free ((void *) acoef);
+ free ((void *) coef);
+ free ((void *) lomega);
+
+
+ oamp = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
+ return oamp;
+}
+
+/*
+
+
+
+compute net rotation from velocities given at r, theta, phi as vel_theta and vel_phi
+
+the routines below are based on code originally from Peter Bird, see
+copyright notice below
+
+
+
+the routine will only properly work for global models if the sampling
+is roughly equal area!
+
+mode: 0 initialize
+ 1 sum
+ 2 solve
+
+*/
+
+
+double determine_netr_tp(float r,float theta,float phi,
+ float velt,float velp,int mode,
+ double *c9,double *omega)
+{
+ float coslat,coslon,sinlat,sinlon,rx,ry,rz,rate,rzu,a,b,c,d,e,f;
+ int i,j,ind[3];
+ double amp,coef[3][3];
+ switch(mode){
+ case 0: /* initialize */
+ for(i=0;i < 9;i++)
+ c9[i] = 0.0;
+ amp = 0.0;
+ break;
+ case 1: /* add this velocity */
+ if((fabs(theta) > 1e-6) &&(fabs(theta-M_PI)>1e-6)){
+ coslat=sin(theta);
+ coslon=cos(phi);
+ sinlat=cos(theta);
+ sinlon=sin(phi);
+
+ rx=coslat*coslon*r;
+ ry=coslat*sinlon*r;
+ rz=sinlat*r;
+
+ rzu=sinlat;
+
+ a = -rz*rzu*sinlon-ry*coslat;
+ b = -rz*coslon;
+ c = rz*rzu*coslon+rx*coslat;
+ d = -rz*sinlon;
+ e = -ry*rzu*coslon+rx*rzu*sinlon;
+ f = ry*sinlon+rx*coslon;
+
+ c9[0] += a*a+b*b;
+ c9[1] += a*c+b*d;
+ c9[2] += a*e+b*f;
+ c9[3] += c*c+d*d;
+ c9[4] += c*e+d*f;
+ c9[5] += e*e+f*f;
+
+ c9[6] += a*velt+b*velp;
+ c9[7] += c*velt+d*velp;
+ c9[8] += e*velt+f*velp;
+ }
+ amp = 0;
+ break;
+ case 2: /* solve */
+ coef[0][0] = c9[0]; /* assemble matrix */
+ coef[0][1] = c9[1];
+ coef[0][2] = c9[2];
+ coef[1][1] = c9[3];
+ coef[1][2] = c9[4];
+ coef[2][2] = c9[5];
+ coef[1][0]=coef[0][1]; /* symmetric */
+ coef[2][0]=coef[0][2];
+ coef[2][1]=coef[1][2];
+ /* */
+ omega[0] = c9[6];
+ omega[1] = c9[7];
+ omega[2] = c9[8];
+ /* solve solution*/
+ hc_ludcmp_3x3(coef,ind);
+ hc_lubksb_3x3(coef,ind,omega);
+ amp = sqrt(omega[0]*omega[0] + omega[1]*omega[1] + omega[2]*omega[2]);
+ break;
+ default:
+ fprintf(stderr,"determine_netr_tp: mode %i undefined\n",mode);
+ parallel_process_termination();
+ break;
+ }
+ return amp;
+}
+
+//
+// subtract a net rotation component from a velocity
+// field given as v_theta (velt) and v_phi (velp) on
+//
+
+void sub_netr(float r,float theta,float phi,float *velt,float *velp, double *omega)
+{
+
+ float coslat,coslon,sinlon,sinlat,rx,ry,rz;
+ float vx,vy,vz,tx,ty,tz,vtheta,pc,px,py,vphi;
+
+ coslat=sin(theta);
+ coslon=cos(phi);
+ sinlat=cos(theta);
+ sinlon=sin(phi);
+
+ rx = coslat*coslon*r;
+ ry = coslat*sinlon*r;
+ rz = sinlat*r;
+
+
+ vx = omega[1]*rz - omega[2]*ry;
+ vy = omega[2]*rx - omega[0]*rz;
+ vz = omega[0]*ry - omega[1]*rx;
+
+ tx = sinlat*coslon; /* theta basis vectors */
+ ty = sinlat*sinlon;
+ tz = -coslat;
+
+ vtheta = vx*tx + vy*ty + vz*tz;
+
+ px = -sinlon; /* phi basis vectors */
+ py = coslon;
+
+ vphi = vx * px + vy * py;
+
+ /* remove */
+ *velt -= vtheta;
+ *velp -= vphi;
+}
+
+
+
+//
+// PROGRAM -OrbScore-: COMPARES OUTPUT FROM -SHELLS-
+// WITH DATA FROM GEODETI// NETWORKS,
+// STRESS DIRECTIONS, FAULT SLIP RATES,
+// SEAFLOOR SPREADING RATES, AND SEISMICITY,
+// AND REPORTS SUMMARY SCALAR SCORES.
+//
+//=========== PART OF THE "SHELLS" PACKAGE OF PROGRAMS===========
+//
+// GIVEN A FINITE ELEMENT GRID FILE, IN THE FORMAT PRODUCED BY
+// -OrbWeave- AND RENUMBERED BY -OrbNumbr-, WITH NODAL DATA
+// ADDED BY -OrbData-, AND NODE-VELOCITY OUTPUT FROM -SHELLS-,
+// COMPUTES A VARIETY OF SCORES OF THE RESULTS.
+//
+// NOTE: Does not contain VISCOS or DIAMND, hence independent
+// of changes made in May 1998, and equally compatible
+// with Old_SHELLS or with improved SHELLS.
+//
+// by
+// Peter Bird
+// Department of Earth and Spcae Sciences,
+// University of California, Los Angeles, California 90095-1567
+// (C) Copyright 1994, 1998, 1999, 2000
+// by Peter Bird and the Regents of
+// the University of California.
+// (For version data see FORMAT 1 below)
+//
+// THIS PROGRAM WAS DEVELOPED WITH SUPPORT FROM THE UNIVERSITY OF
+// CALIFORNIA, THE UNITED STATES GEOLOGI// SURVEY, THE NATIONAL
+// SCIENCE FOUNDATION, AND THE NATIONAL AERONAUTICS AND SPACE
+// ADMINISTRATION.
+// IT IS FREEWARE, AND MAY BE COPIED AND USED WITHOUT CHARGE.
+// IT MAY NOT BE MODIFIED IN A WAY WHICH HIDES ITS ORIGIN
+// OR REMOVES THIS MESSAGE OR THE COPYRIGHT MESSAGE.
+// IT MAY NOT BE RESOLD FOR MORE THAN THE COST OF REPRODUCTION
+// AND MAILING.
+//
+
+
+
+/*
+
+matrix solvers from numerical recipes
+
+ */
+#define NR_TINY 1.0e-20;
+
+void hc_ludcmp_3x3(double a[3][3],int *indx)
+{
+ int i,imax=0,j,k;
+ double big,dum,sum,temp;
+ double vv[3];
+
+ for (i=0;i < 3;i++) {
+ big=0.0;
+ for (j=0;j < 3;j++)
+ if ((temp = fabs(a[i][j])) > big)
+ big=temp;
+ if (fabs(big) < 5e-15) {
+ fprintf(stderr,"hc_ludcmp_3x3: singular matrix in routine, big: %g\n",
+ big);
+ //hc_print_3x3(a,stderr);
+ for(j=0;j<3;j++)
+ fprintf(stderr,"%g %g %g\n",a[j][0],a[j][1],a[j][2]);
+ parallel_process_termination();
+ }
+ vv[i]=1.0/big;
+ }
+ for (j=0;j < 3;j++) {
+ for (i=0;i < j;i++) {
+ sum = a[i][j];
+ for (k=0;k < i;k++)
+ sum -= a[i][k] * a[k][j];
+ a[i][j]=sum;
+ }
+ big=0.0;
+ for (i=j;i < 3;i++) {
+ sum=a[i][j];
+ for (k=0;k < j;k++)
+ sum -= a[i][k] * a[k][j];
+ a[i][j]=sum;
+ if ( (dum = vv[i]*fabs(sum)) >= big) {
+ big=dum;
+ imax=i;
+ }
+ }
+ if (j != imax) {
+ for (k=0;k < 3;k++) {
+ dum = a[imax][k];
+ a[imax][k]=a[j][k];
+ a[j][k]=dum;
+ }
+ vv[imax]=vv[j];
+ }
+ indx[j]=imax;
+ if (fabs(a[j][j]) < 5e-15)
+ a[j][j] = NR_TINY;
+ if (j != 2) {
+ dum=1.0/(a[j][j]);
+ for (i=j+1;i < 3;i++)
+ a[i][j] *= dum;
+ }
+ }
+}
+#undef NR_TINY
+void hc_lubksb_3x3(double a[3][3], int *indx, double *b)
+{
+ int i,ii=0,ip,j;
+ double sum;
+ for (i=0;i < 3;i++) {
+ ip = indx[i];
+ sum = b[ip];
+ b[ip]=b[i];
+ if (ii)
+ for (j=ii-1;j <= i-1;j++)
+ sum -= a[i][j]*b[j];
+ else if (fabs(sum) > 5e-15)
+ ii = i+1;
+ b[i]=sum;
+ }
+ for (i=2;i>=0;i--) {
+ sum=b[i];
+ for (j=i+1;j < 3;j++)
+ sum -= a[i][j]*b[j];
+ b[i] = sum/a[i][i];
+ }
+}
+
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-22 01:00:12 UTC (rev 7864)
@@ -140,7 +140,7 @@
set_sphere_harmonics (E);
- if(E->control.tracer==1) {
+ if(E->control.tracer) {
tracer_initial_settings(E);
(E->problem_tracer_setup)(E);
}
@@ -1316,7 +1316,6 @@
/* delete the geo files */
get_vtk_filename(files,1,E,0);
sprintf(message,"rm %s",files);system(message);
- fprintf(stderr,"%s",message);
if(E->parallel.me == 0){
/* close the log */
fclose(E->output.gzdir.vtk_fp);
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-22 01:00:12 UTC (rev 7864)
@@ -69,6 +69,7 @@
Construct_arrays.c \
Convection.c \
convection_variables.h \
+ Determine_net_rotation.c \
Drive_solvers.c \
drive_solvers.h \
Element_calculations.c \
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-08-22 01:00:12 UTC (rev 7864)
@@ -75,6 +75,7 @@
*/
input_int("gzdir_vtkio",&(E->output.gzdir.vtk_io),"0",m);
+ input_boolean("gzdir_rnr",&(E->output.gzdir.rnr),"off",m); /* remove net rotation? */
E->output.gzdir.vtk_base_init = 0;
E->output.gzdir.vtk_base_save = 1; /* should we save the basis vectors? (memory!) */
//fprintf(stderr,"gzdir: vtkio: %i save basis vectors: %i\n",
Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c 2007-08-22 01:00:12 UTC (rev 7864)
@@ -99,6 +99,10 @@
void gzdir_output_pressure(struct All_variables *, int);
+void sub_netr(float, float, float, float *, float *, double *);
+double determine_model_net_rotation(struct All_variables *,double *);
+
+
void restart_tic_from_gzdir_file(struct All_variables *);
void calc_cbase_at_tp(float , float , float *);
@@ -430,7 +434,8 @@
{
int i, j, k,os;
char output_file[255],output_file2[255],message[255],geo_file[255];
- float cvec[3];
+ float cvec[3],vcorr[3];
+ double omega[3],oamp;
gzFile *gzout;
FILE *fp1;
/* for dealing with several processors */
@@ -455,6 +460,13 @@
E->output.gzdir.vtk_base_init = 1;
}
}
+
+ if(E->output.gzdir.rnr){ /* remove the whole model net rotation */
+ oamp = determine_model_net_rotation(E,omega);
+ if(E->parallel.me == 0)
+ fprintf(stderr,"gzdir_output_velo_temp: removing net rotation: |%8.3e, %8.3e, %8.3e| = %8.3e\n",
+ omega[0],omega[1],omega[2],oamp);
+ }
if((E->output.gzdir.vtk_io == 2) || (E->output.gzdir.vtk_io == 3)){
/*
@@ -524,10 +536,23 @@
fp1 = output_open_mode(output_file,"a");
}
for(k=0,j=1;j <= E->sphere.caps_per_proc;j++,k += os) {
- for(i=1;i<=E->lmesh.nno;i++,k += 9) {
- convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],E->sphere.cap[j].V[1][i],E->sphere.cap[j].V[2][i],
- (E->output.gzdir.vtk_base+k),cvec);
- if(be_write_float_to_file(cvec,3,fp1)!=3)BE_WERROR;
+ if(E->output.gzdir.rnr){
+ /* remove NR */
+ for(i=1;i<=E->lmesh.nno;i++,k += 9) {
+ vcorr[0] = E->sphere.cap[j].V[1][i];
+ vcorr[1] = E->sphere.cap[j].V[2][i];
+ sub_netr(E->sx[j][3][i],E->sx[j][1][i],E->sx[j][2][i],(vcorr+0),(vcorr+1),omega);
+ convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],vcorr[0],vcorr[1],
+ (E->output.gzdir.vtk_base+k),cvec);
+ if(be_write_float_to_file(cvec,3,fp1)!=3)BE_WERROR;
+ }
+ }else{
+ /* regular output */
+ for(i=1;i<=E->lmesh.nno;i++,k += 9) {
+ convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],E->sphere.cap[j].V[1][i],E->sphere.cap[j].V[2][i],
+ (E->output.gzdir.vtk_base+k),cvec);
+ if(be_write_float_to_file(cvec,3,fp1)!=3)BE_WERROR;
+ }
}
}
fclose(fp1);fflush(fp1); /* close file and flush buffer */
@@ -564,7 +589,7 @@
snprintf(output_file2,255,"%s/%d/t.%d.%d",
E->control.data_dir,
cycles,E->parallel.me,cycles);
- }else{ /* vel + T */
+ }else{ /* vel + T for old output */
snprintf(output_file2,255,"%s/%d/velo.%d.%d",
E->control.data_dir,cycles,
E->parallel.me,cycles);
@@ -580,13 +605,26 @@
/* VTK */
for(i=1;i<=E->lmesh.nno;i++)
gzprintf(gzout,"%.6e\n",E->T[j][i]);
- } else {
- /* old */
- for(i=1;i<=E->lmesh.nno;i++)
- gzprintf(gzout,"%.6e %.6e %.6e %.6e\n",
- E->sphere.cap[j].V[1][i],
- E->sphere.cap[j].V[2][i],
- E->sphere.cap[j].V[3][i],E->T[j][i]);
+ } else {
+ /* old velo + T output */
+ if(E->output.gzdir.rnr){
+ /* remove NR */
+ for(i=1;i<=E->lmesh.nno;i++){
+ vcorr[0] = E->sphere.cap[j].V[1][i]; /* vt */
+ vcorr[1] = E->sphere.cap[j].V[2][i]; /* vphi */
+ sub_netr(E->sx[j][3][i],E->sx[j][1][i],E->sx[j][2][i],(vcorr+0),(vcorr+1),omega);
+ gzprintf(gzout,"%.6e %.6e %.6e %.6e\n",
+ vcorr[0],vcorr[1],
+ E->sphere.cap[j].V[3][i],E->T[j][i]);
+
+ }
+ }else{
+ for(i=1;i<=E->lmesh.nno;i++)
+ gzprintf(gzout,"%.6e %.6e %.6e %.6e\n",
+ E->sphere.cap[j].V[1][i],
+ E->sphere.cap[j].V[2][i],
+ E->sphere.cap[j].V[3][i],E->T[j][i]);
+ }
}
}
gzclose(gzout);
@@ -598,15 +636,28 @@
E->control.data_dir,cycles,E->parallel.me,cycles);
gzout = gzdir_output_open(output_file,"w");
for(k=0,j=1;j <= E->sphere.caps_per_proc;j++,k += os) {
- for(i=1;i<=E->lmesh.nno;i++,k += 9) {
- /* convert r,theta,phi vector to x,y,z at base location */
- convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],
- E->sphere.cap[j].V[1][i],
- E->sphere.cap[j].V[2][i],
- (E->output.gzdir.vtk_base+k),cvec);
- /* output of cartesian vector */
- gzprintf(gzout,"%10.4e %10.4e %10.4e\n",
- cvec[0],cvec[1],cvec[2]);
+ if(E->output.gzdir.rnr){
+ /* remove NR */
+ for(i=1;i<=E->lmesh.nno;i++,k += 9) {
+ vcorr[0] = E->sphere.cap[j].V[1][i];
+ vcorr[1] = E->sphere.cap[j].V[2][i];
+ sub_netr(E->sx[j][3][i],E->sx[j][1][i],E->sx[j][2][i],(vcorr+0),(vcorr+1),omega);
+ convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],vcorr[0],vcorr[1],
+ (E->output.gzdir.vtk_base+k),cvec);
+ gzprintf(gzout,"%10.4e %10.4e %10.4e\n",cvec[0],cvec[1],cvec[2]);
+ }
+ }else{
+ /* regular output */
+ for(i=1;i<=E->lmesh.nno;i++,k += 9) {
+ /* convert r,theta,phi vector to x,y,z at base location */
+ convert_pvec_to_cvec(E->sphere.cap[j].V[3][i],
+ E->sphere.cap[j].V[1][i],
+ E->sphere.cap[j].V[2][i],
+ (E->output.gzdir.vtk_base+k),cvec);
+ /* output of cartesian vector */
+ gzprintf(gzout,"%10.4e %10.4e %10.4e\n",
+ cvec[0],cvec[1],cvec[2]);
+ }
}
}
gzclose(gzout);
Modified: mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-22 01:00:12 UTC (rev 7864)
@@ -58,9 +58,9 @@
void myerror(struct All_variables *,char *);
void get_r_spacing_fine(double *, double , double ,int ,double , double ,double , double,
struct All_variables *);
+void xyz2rtp(float ,float ,float ,float *);
-
int get_process_identifier()
{
int pid;
@@ -361,6 +361,15 @@
xout[1] = rst * sin(phi); /* y */
xout[2] = r * cos(theta);
}
+void xyz2rtp(float x,float y,float z,float *rout)
+{
+ float tmp1,tmp2;
+ tmp1 = x*x + y*y;
+ tmp2 = tmp1 + z*z;
+ rout[0] = sqrt(tmp2); /* r */
+ rout[1] = atan2(sqrt(tmp1),z); /* theta */
+ rout[2] = atan2(y,x); /* phi */
+}
/* compute base vectors for conversion of polar to cartesian vectors
Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c 2007-08-22 01:00:12 UTC (rev 7864)
@@ -90,16 +90,21 @@
input_boolean("tracer_enriched",
&(E->control.tracer_enriched),"off",m);
if(E->control.tracer_enriched){
- if(E->composition.ncomp != 1)
- myerror(E,"enriched tracers cannot deal with more than one composition");
if(!E->control.tracer) /* check here so that we can get away
with only one if statement in
Advection_diffusion */
myerror(E,"need to switch on tracers for tracer_enriched");
input_float("Q0_enriched",&(E->control.Q0ER),"0.0",m);
- snprintf(message,100,"using compositionally enriched heating: C = 0: %g C = 1: %g",
+ snprintf(message,100,"using compositionally enriched heating: C = 0: %g C = 1: %g (only one composition!)",
E->control.Q0,E->control.Q0ER);
report(E,message);
+ //
+ // this check doesn't work at this point in the code, and we didn't want to put it into every call to
+ // Advection diffusion
+ //
+ //if(E->composition.ncomp != 1)
+ //myerror(E,"enriched tracers cannot deal with more than one composition yet");
+
}
if(E->control.tracer) {
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-21 22:42:18 UTC (rev 7863)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-22 01:00:12 UTC (rev 7864)
@@ -575,6 +575,8 @@
float *vtk_base;
FILE *vtk_fp;
+ int rnr; /* remove net rotation? */
+
};
struct Output {
More information about the cig-commits
mailing list