[cig-commits] r5838 - in mc/3D/CitcomS/branches/compressible:
CitcomS/Components lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Jan 19 12:14:42 PST 2007
Author: tan2
Date: 2007-01-19 12:14:41 -0800 (Fri, 19 Jan 2007)
New Revision: 5838
Added:
mc/3D/CitcomS/branches/compressible/lib/Regional_tracer_advection.c
mc/3D/CitcomS/branches/compressible/lib/Tracer_setup.c
Removed:
mc/3D/CitcomS/branches/compressible/lib/Tracer_advection.c
Modified:
mc/3D/CitcomS/branches/compressible/CitcomS/Components/IC.py
mc/3D/CitcomS/branches/compressible/CitcomS/Components/Visc.py
mc/3D/CitcomS/branches/compressible/lib/Full_read_input_from_files.c
mc/3D/CitcomS/branches/compressible/lib/Instructions.c
mc/3D/CitcomS/branches/compressible/lib/Makefile.am
mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/branches/compressible/lib/Regional_read_input_from_files.c
mc/3D/CitcomS/branches/compressible/lib/Regional_version_dependent.c
mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c
mc/3D/CitcomS/branches/compressible/lib/tracer_defs.h
mc/3D/CitcomS/branches/compressible/lib/viscosity_descriptions.h
mc/3D/CitcomS/branches/compressible/module/misc.c
mc/3D/CitcomS/branches/compressible/module/setProperties.c
Log:
Porting r5744-5837 from trunk
Modified: mc/3D/CitcomS/branches/compressible/CitcomS/Components/IC.py
===================================================================
--- mc/3D/CitcomS/branches/compressible/CitcomS/Components/IC.py 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/CitcomS/Components/IC.py 2007-01-19 20:14:41 UTC (rev 5838)
@@ -116,9 +116,11 @@
perturblayer = pyre.inventory.slice("perturblayer", default=[5])
perturbmag = pyre.inventory.list("perturbmag", default=[0.05])
- # for tic_method=2
+ # for tic_method=1 or 2
half_space_age = pyre.inventory.float("half_space_age", default=40,
validator=pyre.inventory.greater(1e-3))
+
+ # for tic_method=2
blob_center = pyre.inventory.list("blob_center", default=[-999., -999., -999.])
blob_radius = pyre.inventory.float("blob_radius", default=0.063)
blob_dT = pyre.inventory.float("blob_dT", default=0.18)
Modified: mc/3D/CitcomS/branches/compressible/CitcomS/Components/Visc.py
===================================================================
--- mc/3D/CitcomS/branches/compressible/CitcomS/Components/Visc.py 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/CitcomS/Components/Visc.py 2007-01-19 20:14:41 UTC (rev 5838)
@@ -84,6 +84,14 @@
sdepv_expt = pyre.inventory.list("sdepv_expt", default=[1, 1, 1, 1])
sdepv_misfit = pyre.inventory.float("sdepv_misfit", default=0.02)
+ low_visc_channel = pyre.inventory.bool("low_visc_channel", default=False)
+ low_visc_wedge = pyre.inventory.bool("low_visc_wedge", default=False)
+
+ lv_min_radius = pyre.inventory.float("lv_min_radius", default=0.9764)
+ lv_max_radius = pyre.inventory.float("lv_max_radius", default=0.9921)
+ lv_channel_thickness = pyre.inventory.float("lv_channel_thickness", default=0.0047)
+ lv_reduction = pyre.inventory.float("lv_reduction", default=0.5)
+
VMIN = pyre.inventory.bool("VMIN", default=False)
visc_min = pyre.inventory.float("visc_min", default=1.0e-3)
Modified: mc/3D/CitcomS/branches/compressible/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Full_read_input_from_files.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Full_read_input_from_files.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -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>
@@ -166,8 +166,8 @@
E->sphere.cap[m].VB[3][nodel] = 0.0;
}
else { /* negative ages - don't do the interpolation */
- E->sphere.cap[m].VB[1][nodel] = VB1[1][nodeg];
- E->sphere.cap[m].VB[2][nodel] = VB1[2][nodeg];
+ E->sphere.cap[m].VB[1][nodel] = VB1[1][nodeg] * E->data.scalev;
+ E->sphere.cap[m].VB[2][nodel] = VB1[2][nodeg] * E->data.scalev;
E->sphere.cap[m].VB[3][nodel] = 0.0;
}
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -254,7 +254,7 @@
set_mg_defaults(E);}
else
{ if (E->parallel.me==0) fprintf(stderr,"Unable to determine how to solve, specify Solver=VALID_OPTION \n");
- exit(0);
+ parallel_process_termination();
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Makefile.am 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Makefile.am 2007-01-19 20:14:41 UTC (rev 5838)
@@ -106,8 +106,8 @@
Sphere_util.c \
Stokes_flow_Incomp.c \
Topo_gravity.c \
- Tracer_advection.c \
tracer_defs.h \
+ Tracer_setup.c \
viscosity_descriptions.h \
Viscosity_structures.c \
Full_boundary_conditions.c \
@@ -125,6 +125,7 @@
Regional_read_input_from_files.c \
Regional_solver.c \
Regional_sphere_related.c \
+ Regional_tracer_advection.c \
Regional_version_dependent.c
EXTRA_DIST = \
Modified: mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -184,7 +184,7 @@
ptr = malloc((size_t)bytes);
if (ptr == (void *)NULL) {
fprintf(stderr,"Memory: cannot allocate another %d bytes \n(line %d of file %s)\n",bytes,line,file);
- exit(0);
+ parallel_process_termination();
}
return(ptr);
@@ -246,7 +246,7 @@
fprintf(E->fp,"Unable to open the required file `%s' (this is fatal)",filename);
fflush(E->fp);
- exit(1);
+ parallel_process_termination();
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Regional_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Regional_read_input_from_files.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Regional_read_input_from_files.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -136,8 +136,8 @@
E->sphere.cap[1].VB[3][nodel] = 0.0;
}
else { /* negative ages - don't do the interpolation */
- E->sphere.cap[1].VB[1][nodel] = VB1[1][nodeg];
- E->sphere.cap[1].VB[2][nodel] = VB1[2][nodeg];
+ E->sphere.cap[1].VB[1][nodel] = VB1[1][nodeg]*E->data.scalev;
+ E->sphere.cap[1].VB[2][nodel] = VB1[2][nodeg]*E->data.scalev;
E->sphere.cap[1].VB[3][nodel] = 0.0;
}
}
Added: mc/3D/CitcomS/branches/compressible/lib/Regional_tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Regional_tracer_advection.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Regional_tracer_advection.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -0,0 +1,822 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<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>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+/*
+
+Regional_tracer_advection.c
+
+A program which initiates the distribution of tracers
+and advects those tracers in a time evolving velocity field.
+Called and used from the CitCOM finite element code.
+Written 2/96 M. Gurnis for Citcom in cartesian geometry
+Modified by Lijie in 1998 and by Vlad and Eh in 2005 for CitcomS
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <mpi.h>
+#include <math.h>
+#include <sys/types.h>
+#ifdef HAVE_MALLOC_H
+#include <malloc.h>
+#endif
+
+#include "element_definitions.h"
+#include "global_defs.h"
+
+static void mean_elem_coord(struct All_variables *E);
+
+void regional_tracer_setup(struct All_variables *E)
+{
+ int i,j,k,node;
+ int m,ntr;
+ int n_x,n_y,n_z;
+ int node1,node2,node3,node4,node5,node6,node7,node8;
+ int local_element;
+ float THETA_LOC_ELEM,FI_LOC_ELEM,R_LOC_ELEM;
+ float idummy,xdummy,ydummy,zdummy;
+ FILE *fp;
+ int nox,noy,noz;
+ char output_file[255];
+ MPI_Comm world;
+ MPI_Status status;
+
+ n_x=0;
+ n_y=0;
+ n_z=0;
+
+ nox=E->lmesh.nox;
+ noy=E->lmesh.noy;
+ noz=E->lmesh.noz;
+
+ sprintf(output_file,"%s",E->control.tracer_file);
+ fp=fopen(output_file,"r");
+ if (fp == NULL) {
+ fprintf(E->fp,"(Tracer_advection #1) Cannot open %s\n", output_file);
+ exit(8);
+ }
+ fscanf(fp,"%d",&(E->Tracer.NUM_TRACERS));
+
+ E->Tracer.tracer_x=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ E->Tracer.tracer_y=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ E->Tracer.tracer_z=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ E->Tracer.itcolor=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+
+ E->Tracer.x_space=(float*) malloc((nox+1)*sizeof(float));
+ E->Tracer.y_space=(float*) malloc((noy+1)*sizeof(float));
+ E->Tracer.z_space=(float*) malloc((noz+1)*sizeof(float));
+
+ E->Tracer.LOCAL_ELEMENT=(int*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(int));
+
+ /* for rheology stuff */
+
+ E->Tracer.THETA_LOC_ELEM=(float*) malloc((E->lmesh.nno+1)*sizeof(float));
+ E->Tracer.FI_LOC_ELEM=(float*) malloc((E->lmesh.nno+1)*sizeof(float));
+ E->Tracer.R_LOC_ELEM=(float*) malloc((E->lmesh.nno+1)*sizeof(float));
+
+ E->Tracer.THETA_LOC_ELEM_T=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ E->Tracer.FI_LOC_ELEM_T=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ E->Tracer.R_LOC_ELEM_T=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+
+
+
+ /***comment by Vlad 03/15/2005
+ each processor holds its own number of tracers
+ ***/
+
+ ntr=1;
+ for(i=1;i<=E->Tracer.NUM_TRACERS;i++) {
+ fscanf(fp,"%f %f %f %f", &idummy, &xdummy, &ydummy, &zdummy);
+ if(xdummy >= E->sx[1][1][1] && xdummy <= E->sx[1][1][nox*noy*noz]) {
+ if(ydummy >= E->sx[1][2][1] && ydummy <= E->sx[1][2][nox*noy*noz]) {
+ if(zdummy >= E->sx[1][3][1] && zdummy <= E->sx[1][3][nox*noy*noz]) {
+ E->Tracer.itcolor[ntr]=idummy;
+ E->Tracer.tracer_x[ntr]=xdummy;
+ E->Tracer.tracer_y[ntr]=ydummy;
+ E->Tracer.tracer_z[ntr]=zdummy;
+ ntr++;
+ }
+ }
+ }
+ }
+
+ /***comment by Vlad 3/30/2005
+ E->Tracer.LOCAL_NUM_TRACERS is the initial number
+ of tracers in each processor
+ ***/
+
+ E->Tracer.LOCAL_NUM_TRACERS=ntr-1;
+
+
+
+ /***comment by Vlad 1/26/2005
+ reading the local mesh coordinate
+ ***/
+
+
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(i=1;i<=nox;i++)
+ {
+ j=1;
+ k=1;
+ node=k+(i-1)*noz+(j-1)*nox*noz;
+ E->Tracer.x_space[i]=E->sx[m][1][node];
+ }
+
+ for(j=1;j<=noy;j++)
+ {
+ i=1;
+ k=1;
+ node=k+(i-1)*noz+(j-1)*nox*noz;
+ E->Tracer.y_space[j]=E->sx[m][2][node];
+ }
+
+ for(k=1;k<=noz;k++)
+ {
+ i=1;
+ j=1;
+ node=k+(i-1)*noz+(j-1)*nox*noz;
+ E->Tracer.z_space[k]=E->sx[m][3][node];
+ }
+
+ } /* end of m */
+
+ /***comment by Vlad 04/20/2006
+ reading the local element eight nodes coordinate,
+ then computing the mean element coordinates
+ ***/
+
+
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(j=1;j<=(noy-1);j++) {
+ for(i=1;i<=(nox-1);i++) {
+ for(k=1;k<=(noz-1);k++) {
+ n_x=i;
+ n_y=j;
+ n_z=k;
+
+ node1 = n_z + (n_x-1)*noz + (n_y-1)*noz*nox;
+ node2 = n_z + n_x*noz + (n_y-1)*noz*nox;
+ node3 = n_z+1 + (n_x-1)*noz + (n_y-1)*noz*nox;
+ node4 = n_z+1 + n_x*noz + (n_y-1)*noz*nox;
+ node5 = n_z + (n_x-1)*noz + n_y*noz*nox;
+ node6 = n_z + n_x*noz + n_y*noz*nox;
+ node7 = n_z+1 + (n_x-1)*noz + n_y*noz*nox;
+ node8 = n_z+1 + n_x*noz + n_y*noz*nox;
+
+ /* for rheology stuff */
+ local_element=node1-(n_x-1)-(n_y-1)*(nox+noz-1);
+
+ E->Tracer.THETA_LOC_ELEM[local_element]=(E->sx[m][1][node1]+E->sx[m][1][node2])/2;
+ E->Tracer.FI_LOC_ELEM[local_element]=(E->sx[m][2][node1]+E->sx[m][2][node5])/2;
+ E->Tracer.R_LOC_ELEM[local_element]=(E->sx[m][3][node1]+E->sx[m][3][node3])/2;
+
+ }
+ }
+ }
+ } /* end of m */
+
+ mean_elem_coord(E);
+
+ return;
+
+}
+
+
+
+
+void regional_tracer_advection(struct All_variables *E)
+{
+ int i,j,k,l,m,n,o,p;
+ int n_x,n_y,n_z;
+ int node1,node2,node3,node4,node5,node6,node7,node8;
+ int nno,nox,noy,noz;
+ int iteration;
+ float x_tmp, y_tmp, z_tmp;
+ float volume, tr_dx, tr_dy, tr_dz, dx, dy, dz;
+ float w1,w2,w3,w4,w5,w6,w7,w8;
+ float tr_v[NCS][4];
+ MPI_Comm world;
+ MPI_Status status[4];
+ MPI_Status status_count;
+ MPI_Request request[4];
+ float xmin,xmax,ymin,ymax,zmin,zmax;
+
+ float x_global_min,x_global_max,y_global_min,y_global_max,z_global_min,z_global_max;
+
+
+ float *tr_color_1,*tr_x_1,*tr_y_1,*tr_z_1;
+ float *tr_color_new[13],*tr_x_new[13],*tr_y_new[13],*tr_z_new[13];
+ int *Left_proc_list,*Right_proc_list;
+ int *jump_new,*count_new;
+ int *jj;
+
+ int proc;
+ int Previous_num_tracers,Current_num_tracers;
+ int locx,locy,locz;
+ int left,right,up,down,back,front;
+ int temp_tracers;
+
+
+ world=E->parallel.world;
+
+
+ nox=E->lmesh.nox;
+ noy=E->lmesh.noy;
+ noz=E->lmesh.noz;
+ nno=nox*noy*noz;
+
+ Left_proc_list=(int*) malloc(6*sizeof(int));
+ Right_proc_list=(int*) malloc(6*sizeof(int));
+ jump_new=(int*) malloc(6*sizeof(int));
+ count_new=(int*) malloc(6*sizeof(int));
+ jj=(int*) malloc(6*sizeof(int));
+
+ tr_x_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ tr_y_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ tr_z_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ tr_color_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+
+ for(i=0;i<=11;i++){
+ tr_color_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ tr_x_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ tr_y_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ tr_z_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
+ }
+
+
+ /*** comment by Vlad 3/25/2005
+ This part of code gets the bounding box of the local mesh.
+ ***/
+
+ xmin=E->Tracer.x_space[1];
+ xmax=E->Tracer.x_space[nox];
+ ymin=E->Tracer.y_space[1];
+ ymax=E->Tracer.y_space[noy];
+ zmin=E->Tracer.z_space[1];
+ zmax=E->Tracer.z_space[noz];
+
+ /*fprintf(stderr,"%d %d\n", E->parallel.nprocx, E->parallel.loc2proc_map[0][0][1][0]);*/
+
+ /*** comment by Tan2 1/25/2005
+ Copy the velocity array.
+ ***/
+
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(i=1;i<=nno;i++)
+ for(j=1;j<=3;j++) {
+ E->GV[m][j][i]=E->sphere.cap[m].V[j][i];
+ }
+ }
+
+ /*** comment by vlad 03/17/2005
+ advecting tracers in each processor
+ ***/
+
+
+ for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
+
+ n_x=0;
+ n_y=0;
+ n_z=0;
+
+
+ /* mid point method uses 2 iterations */
+
+ x_tmp=E->Tracer.tracer_x[n];
+ y_tmp=E->Tracer.tracer_y[n];
+ z_tmp=E->Tracer.tracer_z[n];
+
+ for(iteration=1;iteration<=2;iteration++)
+ {
+
+
+ /*** comment by Tan2 1/25/2005
+ Find the element that contains the tracer.
+
+ nodex n_x n_x+1
+ | * |
+ <----------->
+ tr_dx
+
+ <-------------------------------->
+ dx
+ ***/
+
+ for(i=1;i<nox;i++) {
+ if(x_tmp >= E->Tracer.x_space[i] && x_tmp <= E->Tracer.x_space[i+1]) {
+ tr_dx=x_tmp-E->Tracer.x_space[i];
+ dx=E->Tracer.x_space[i+1]-E->Tracer.x_space[i];
+ n_x=i;
+
+ E->Tracer.THETA_LOC_ELEM_T[n]=(E->Tracer.x_space[i+1]+E->Tracer.x_space[i])/2;
+ }
+ }
+
+ for(j=1;j<noy;j++) {
+ if(y_tmp >= E->Tracer.y_space[j] && y_tmp <= E->Tracer.y_space[j+1]) {
+ tr_dy=y_tmp-E->Tracer.y_space[j];
+ dy=E->Tracer.y_space[j+1]-E->Tracer.y_space[j];
+ n_y=j;
+
+ E->Tracer.FI_LOC_ELEM_T[n]=(E->Tracer.y_space[j+1]+E->Tracer.y_space[j])/2;
+
+ }
+ }
+
+ for(k=1;k<noz;k++) {
+ if(z_tmp >= E->Tracer.z_space[k] && z_tmp <= E->Tracer.z_space[k+1]) {
+ tr_dz=z_tmp-E->Tracer.z_space[k];
+ dz=E->Tracer.z_space[k+1]-E->Tracer.z_space[k];
+ n_z=k;
+
+ E->Tracer.R_LOC_ELEM_T[n]=(E->Tracer.z_space[k+1]+E->Tracer.z_space[k])/2;
+
+ }
+ }
+
+ //fprintf(stderr,"tracer: %d %f %f %f\n",n,E->Tracer.THETA_LOC_ELEM_T[n],E->Tracer.FI_LOC_ELEM_T[n],E->Tracer.R_LOC_ELEM_T[n]);
+
+
+ /*** comment by Tan2 1/25/2005
+ Calculate shape functions from tr_dx, tr_dy, tr_dz
+ This assumes linear element
+ ***/
+
+
+ /* compute volumetic weighting functions */
+ w1=tr_dx*tr_dz*tr_dy;
+ w2=(dx-tr_dx)*tr_dz*tr_dy;
+ w3=tr_dx*(dz-tr_dz)*tr_dy;
+ w4=(dx-tr_dx)*(dz-tr_dz)*tr_dy;
+ w5=tr_dx*tr_dz*(dy-tr_dy);
+ w6=(dx-tr_dx)*tr_dz*(dy-tr_dy);
+ w7=tr_dx*(dz-tr_dz)*(dy-tr_dy);
+ w8=(dx-tr_dx)*(dz-tr_dz)*(dy-tr_dy);
+
+
+ volume=dx*dz*dy;
+
+
+ /*** comment by Tan2 1/25/2005
+ Calculate the 8 node numbers of current element
+ ***/
+
+ node1 = n_z + (n_x-1)*noz + (n_y-1)*noz*nox;
+ node2 = n_z + n_x*noz + (n_y-1)*noz*nox;
+ node3 = n_z+1 + (n_x-1)*noz + (n_y-1)*noz*nox;
+ node4 = n_z+1 + n_x*noz + (n_y-1)*noz*nox;
+ node5 = n_z + (n_x-1)*noz + n_y*noz*nox;
+ node6 = n_z + n_x*noz + n_y*noz*nox;
+ node7 = n_z+1 + (n_x-1)*noz + n_y*noz*nox;
+ node8 = n_z+1 + n_x*noz + n_y*noz*nox;
+
+
+ /* printf("%d %d %d %d %d %d %d %d %d\n",E->parallel.me, node1,node2,node3,node4,node5,node6,node7,node8);
+ //printf("%d %f %f %f %f %f %f %f %f\n", E->parallel.me, E->GV[1][2][node1], E->GV[1][2][node2], E->GV[1][2][node3], E->GV[1][2][node4], E->GV[1][2][node5], E->GV[1][2][node6], E->GV[1][2][node7], E->GV[1][2][node8]);
+ */
+
+ /*** comment by Tan2 1/25/2005
+ Interpolate the velocity on the tracer position
+ ***/
+
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(j=1;j<=3;j++) {
+ tr_v[m][j]=w8*E->GV[m][j][node1]
+ +w7*E->GV[m][j][node2]
+ +w6*E->GV[m][j][node3]
+ +w5*E->GV[m][j][node4]
+ +w4*E->GV[m][j][node5]
+ +w3*E->GV[m][j][node6]
+ +w2*E->GV[m][j][node7]
+ +w1*E->GV[m][j][node8];
+ tr_v[m][j]=tr_v[m][j]/volume;
+
+ }
+
+
+
+ E->Tracer.LOCAL_ELEMENT[n]=node1-(n_x-1)-(n_y-1)*(nox+noz-1);
+
+
+
+
+ //fprintf(stderr,"%s %d %s %d %f %f %f %f\n", "The tracer no:", n,"is in element no:", E->Tracer.LOCAL_ELEMENT[n], E->Tracer.y_space[n_y], E->Tracer.tracer_y[n], E->Tracer.FI_LOC_ELEM[515], E->Tracer.y_space[n_y+1]);
+
+
+
+
+
+ /*** comment by Tan2 1/25/2005
+ advect tracer using mid-point method (2nd order accuracy)
+ ***/
+
+ /* mid point method */
+
+ if(iteration == 1) {
+ x_tmp = x_tmp + (E->advection.timestep/2.0)*tr_v[m][1]/E->Tracer.z_space[n_z];
+ y_tmp = y_tmp + (E->advection.timestep/2.0)*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
+ z_tmp = z_tmp + (E->advection.timestep/2.0)*tr_v[m][3];
+ }
+ if( iteration == 2) {
+ E->Tracer.tracer_x[n] += E->advection.timestep*tr_v[m][1]/E->Tracer.z_space[n_z];
+ E->Tracer.tracer_y[n] += E->advection.timestep*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
+ E->Tracer.tracer_z[n] += E->advection.timestep*tr_v[m][3];
+
+
+ //fprintf(stderr,"%d %d %f %f %f %f %f %f\n", E->parallel.me, E->monitor.solution_cycles, E->Tracer.tracer_x[n],E->Tracer.tracer_y[n],E->Tracer.tracer_z[n], tr_v[m][1],tr_v[m][2],tr_v[m][3]);
+
+
+ }
+
+
+
+ } /* end of m */
+
+
+ } /* end of iteration loop */
+
+
+ /*** Comment by Vlad 12/15/2005
+ Put the tracers back in the box if they go out
+ ***/
+
+ /*** Comment by Vlad 12/15/2005
+ get the bounding box of the global mesh
+ ***/
+
+ x_global_min = E->control.theta_min;
+ x_global_max = E->control.theta_max;
+ y_global_min = E->control.fi_min;
+ y_global_max = E->control.fi_max;
+ z_global_min = E->sphere.ri;
+ z_global_max = E->sphere.ro;
+
+ //printf("%f %f %f %f %f %f\n", E->sphere.cap[1].theta[1],E->sphere.cap[1].theta[3],E->sphere.cap[1].fi[1],E->sphere.cap[1].fi[3],E->sphere.ri,E->sphere.ro);
+
+ if(E->Tracer.tracer_x[n] > x_global_max)
+ E->Tracer.tracer_x[n] = x_global_max;
+ if(E->Tracer.tracer_x[n] < x_global_min)
+ E->Tracer.tracer_x[n] = x_global_min;
+ if(E->Tracer.tracer_y[n] > y_global_max)
+ E->Tracer.tracer_y[n] = y_global_max;
+ if(E->Tracer.tracer_y[n] < y_global_min)
+ E->Tracer.tracer_y[n] = y_global_min;
+ if(E->Tracer.tracer_z[n] > z_global_max)
+ E->Tracer.tracer_z[n] = z_global_max;
+ if(E->Tracer.tracer_z[n] < z_global_min)
+ E->Tracer.tracer_z[n] = z_global_min;
+
+
+
+ }/* end of tracer loop */
+
+ /*** Comment by Vlad 3/25/2005
+ MPI for the tracer-advection code
+ ***/
+
+
+ m = 0;
+
+ locx = E->parallel.me_loc[1];
+ locy = E->parallel.me_loc[2];
+ locz = E->parallel.me_loc[3];
+
+ /* Am I the left-most proc.? If not, who is on my left? */
+ if (locy == 0)
+ left = -1;
+ else
+ left = E->parallel.loc2proc_map[m][locx][locy-1][locz];
+
+ /* Am I the right-most proc.? If not, who is on my right? */
+ if (locy == E->parallel.nprocy-1)
+ right = -1;
+ else
+ right = E->parallel.loc2proc_map[m][locx][locy+1][locz];
+
+ /* Am I the lower-most proc.? If not, who is beneath me? */
+ if (locz == 0)
+ down = -1;
+ else
+ down = E->parallel.loc2proc_map[m][locx][locy][locz-1];
+
+ /* Am I the upper-most proc.? If not, who is above me? */
+ if (locz == E->parallel.nprocz-1)
+ up = -1;
+ else
+ up = E->parallel.loc2proc_map[m][locx][locy][locz+1];
+
+ /* Am I the back-most proc.? If not, who is behind me? */
+ if (locx == 0)
+ back = -1;
+ else
+ back = E->parallel.loc2proc_map[m][locx-1][locy][locz];
+
+ /* Am I the front-most proc.? If not, who is in front of me? */
+ if (locx == E->parallel.nprocx-1)
+ front = -1;
+ else
+ front = E->parallel.loc2proc_map[m][locx+1][locy][locz];
+
+
+ Left_proc_list[0]=left;
+ Left_proc_list[1]=right;
+ Left_proc_list[2]=down;
+ Left_proc_list[3]=up;
+ Left_proc_list[4]=back;
+ Left_proc_list[5]=front;
+
+ Right_proc_list[0]=right;
+ Right_proc_list[1]=left;
+ Right_proc_list[2]=up;
+ Right_proc_list[3]=down;
+ Right_proc_list[4]=front;
+ Right_proc_list[5]=back;
+
+ jump_new[0]=0;
+ jump_new[1]=0;
+ jump_new[2]=0;
+ jump_new[3]=0;
+ jump_new[4]=0;
+ jump_new[5]=0;
+
+ count_new[0]=0;
+ count_new[1]=0;
+ count_new[2]=0;
+ count_new[3]=0;
+ count_new[4]=0;
+ count_new[5]=0;
+
+ jj[0]=1;
+ jj[1]=0;
+ jj[2]=3;
+ jj[3]=2;
+ jj[4]=5;
+ jj[5]=4;
+
+ temp_tracers=0;
+ Current_num_tracers=0;
+
+ for(i=0;i<=11;i++){
+ for(j=0;j<=E->Tracer.NUM_TRACERS;j++){
+ tr_color_new[i][j]=999;
+ tr_x_new[i][j]=999;
+ tr_y_new[i][j]=999;
+ tr_z_new[i][j]=999;
+
+ tr_color_1[j]=999;
+ tr_x_1[j]=999;
+ tr_y_1[j]=999;
+ tr_z_1[j]=999;
+ }
+ }
+
+
+ i=0;
+ j=0;
+ k=0;
+ l=0;
+ m=0;
+ o=0;
+ p=0;
+
+
+ for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++){
+
+ if(E->Tracer.tracer_y[n]>ymax) {
+ /* excluding Nan */
+ if(E->Tracer.tracer_y[n]+100 != 100) {
+ tr_color_new[0][i]=E->Tracer.itcolor[n];
+ tr_x_new[0][i]=E->Tracer.tracer_x[n];
+ tr_y_new[0][i]=E->Tracer.tracer_y[n];
+ tr_z_new[0][i]=E->Tracer.tracer_z[n];
+ i++;
+ jump_new[0]=i;
+ }
+ }
+ else if(E->Tracer.tracer_y[n]<ymin) {
+ if(E->Tracer.tracer_y[n]+100 != 100) {
+ tr_color_new[1][j]=E->Tracer.itcolor[n];
+ tr_x_new[1][j]=E->Tracer.tracer_x[n];
+ tr_y_new[1][j]=E->Tracer.tracer_y[n];
+ tr_z_new[1][j]=E->Tracer.tracer_z[n];
+ j++;
+ jump_new[1]=j;
+ }
+ }
+ else if(E->Tracer.tracer_z[n]>zmax) {
+ if(E->Tracer.tracer_z[n]+100 != 100) {
+ tr_color_new[2][k]=E->Tracer.itcolor[n];
+ tr_x_new[2][k]=E->Tracer.tracer_x[n];
+ tr_y_new[2][k]=E->Tracer.tracer_y[n];
+ tr_z_new[2][k]=E->Tracer.tracer_z[n];
+ k++;
+ jump_new[2]=k;
+ }
+ }
+ else if(E->Tracer.tracer_z[n]<zmin) {
+ if(E->Tracer.tracer_z[n]+100 != 100) {
+ tr_color_new[3][l]=E->Tracer.itcolor[n];
+ tr_x_new[3][l]=E->Tracer.tracer_x[n];
+ tr_y_new[3][l]=E->Tracer.tracer_y[n];
+ tr_z_new[3][l]=E->Tracer.tracer_z[n];
+ l++;
+ jump_new[3]=l;
+ }
+ }
+
+ else if(E->Tracer.tracer_x[n]>xmax) {
+ if(E->Tracer.tracer_x[n]+100 != 100) {
+ tr_color_new[4][m]=E->Tracer.itcolor[n];
+ tr_x_new[4][m]=E->Tracer.tracer_x[n];
+ tr_y_new[4][m]=E->Tracer.tracer_y[n];
+ tr_z_new[4][m]=E->Tracer.tracer_z[n];
+ m++;
+ jump_new[4]=m;
+ }
+ }
+ else if(E->Tracer.tracer_x[n]<xmin) {
+ if(E->Tracer.tracer_x[n]+100 != 100) {
+ tr_color_new[5][o]=E->Tracer.itcolor[n];
+ tr_x_new[5][o]=E->Tracer.tracer_x[n];
+ tr_y_new[5][o]=E->Tracer.tracer_y[n];
+ tr_z_new[5][o]=E->Tracer.tracer_z[n];
+ o++;
+ jump_new[5]=o;
+ }
+ }
+
+ else {
+ tr_color_1[p]=E->Tracer.itcolor[n];
+ tr_x_1[p]=E->Tracer.tracer_x[n];
+ tr_y_1[p]=E->Tracer.tracer_y[n];
+ tr_z_1[p]=E->Tracer.tracer_z[n];
+ p++;
+ }
+ }
+
+ Previous_num_tracers=E->Tracer.LOCAL_NUM_TRACERS;
+ Current_num_tracers=Previous_num_tracers-jump_new[0]-jump_new[1]-jump_new[2]-jump_new[3]-jump_new[4]-jump_new[5];
+
+ /* compact the remaining tracer */
+ for(p=1;p<=Current_num_tracers;p++){
+ E->Tracer.itcolor[p]=tr_color_1[p-1];
+ E->Tracer.tracer_x[p]=tr_x_1[p-1];
+ E->Tracer.tracer_y[p]=tr_y_1[p-1];
+ E->Tracer.tracer_z[p]=tr_z_1[p-1];
+ }
+
+
+ for(i=0;i<=5;i++){
+
+ j=jj[i];
+
+ if (Left_proc_list[i] >= 0) {
+ proc=Left_proc_list[i];
+ MPI_Irecv(tr_color_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 11+i, world, &request[0]);
+ MPI_Irecv(tr_x_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 12+i, world, &request[1]);
+ MPI_Irecv(tr_y_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 13+i, world, &request[2]);
+ MPI_Irecv(tr_z_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 14+i, world, &request[3]);
+ }
+
+ if (Right_proc_list[i] >= 0) {
+ proc=Right_proc_list[i];
+ MPI_Send(tr_color_new[i], jump_new[i], MPI_FLOAT, proc, 11+i, world);
+ MPI_Send(tr_x_new[i], jump_new[i], MPI_FLOAT, proc, 12+i, world);
+ MPI_Send(tr_y_new[i], jump_new[i], MPI_FLOAT, proc, 13+i, world);
+ MPI_Send(tr_z_new[i], jump_new[i], MPI_FLOAT, proc, 14+i, world);
+ }
+
+ if (Left_proc_list[i] >= 0) {
+ MPI_Waitall(4, request, status);
+ status_count = status[0];
+ MPI_Get_count(&status_count, MPI_FLOAT, &count_new[i]);
+ }
+
+
+ temp_tracers=temp_tracers+count_new[i]-jump_new[i];
+ E->Tracer.LOCAL_NUM_TRACERS=Previous_num_tracers+temp_tracers;
+
+
+ /* append the tracers */
+
+ if(i <= 1){
+ for(n=Current_num_tracers+count_new[j];n<=Current_num_tracers+count_new[i]+count_new[j]-1;n++) {
+ m=Current_num_tracers+count_new[j];
+ E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
+ E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
+ E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
+ E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
+
+ }
+ }
+
+
+ else if (i <= 3) {
+ for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];n<=Current_num_tracers+count_new[0]+count_new[1]+count_new[i]+count_new[j]-1;n++) {
+ m=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];
+ E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
+ E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
+ E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
+ E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
+
+ }
+ }
+
+ else {
+ for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];n<=E->Tracer.LOCAL_NUM_TRACERS-1;n++) {
+ m=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];
+ E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
+ E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
+ E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
+ E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
+
+ }
+ }
+
+
+ }
+
+
+ free (tr_color_1);
+ free (tr_x_1);
+ free (tr_y_1);
+ free (tr_z_1);
+ for(i=0;i<=11;i++) {
+ free (tr_color_new[i]);
+ free (tr_x_new[i]);
+ free (tr_y_new[i]);
+ free (tr_z_new[i]);
+ }
+
+
+ return;
+}
+
+
+
+/*avoid 1st time step problem with E->Tracer.THETA_LOC_ELEM_T[n],E->Tracer.FI_LOC_ELEM_T[n],E->Tracer.R_LOC_ELEM_T[n] */
+
+static void mean_elem_coord(struct All_variables *E)
+{
+ int n, i, j, k;
+ float x_tmp, y_tmp, z_tmp;
+
+ for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
+ x_tmp=E->Tracer.tracer_x[n];
+ y_tmp=E->Tracer.tracer_y[n];
+ z_tmp=E->Tracer.tracer_z[n];
+
+ for(i=1;i<E->lmesh.nox;i++) {
+ if(x_tmp >= E->Tracer.x_space[i] && x_tmp <= E->Tracer.x_space[i+1]) {
+ E->Tracer.THETA_LOC_ELEM_T[n]=(E->Tracer.x_space[i+1]+E->Tracer.x_space[i])/2;
+ }
+ }
+
+ for(j=1;j<E->lmesh.noy;j++) {
+ if(y_tmp >= E->Tracer.y_space[j] && y_tmp <= E->Tracer.y_space[j+1]) {
+ E->Tracer.FI_LOC_ELEM_T[n]=(E->Tracer.y_space[j+1]+E->Tracer.y_space[j])/2;
+
+ }
+ }
+
+ for(k=1;k<E->lmesh.noz;k++) {
+ if(z_tmp >= E->Tracer.z_space[k] && z_tmp <= E->Tracer.z_space[k+1]) {
+ E->Tracer.R_LOC_ELEM_T[n]=(E->Tracer.z_space[k+1]+E->Tracer.z_space[k])/2;
+
+ }
+ }
+ }
+ return;
+}
Modified: mc/3D/CitcomS/branches/compressible/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Regional_version_dependent.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Regional_version_dependent.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -387,6 +387,16 @@
}
else if (E->convection.tic_method == 1) {
+ /* set up a top thermal boundary layer */
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=noy;i++)
+ for(j=1;j<=nox;j++)
+ for(k=1;k<=noz;k++) {
+ node=k+(j-1)*noz+(i-1)*nox*noz;
+ r1=E->sx[m][3][node];
+ temp = 0.2*(E->sphere.ro-r1) * 0.5/sqrt(E->convection.half_space_age/E->data.scalet);
+ E->T[m][node] = E->control.TBCbotval*erf(temp);
+ }
}
else if (E->convection.tic_method == 2) {
Deleted: mc/3D/CitcomS/branches/compressible/lib/Tracer_advection.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Tracer_advection.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Tracer_advection.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -1,698 +0,0 @@
-/*
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
- *<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>
- *
- *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- */
-/*
-
- Tracer_advection.c
-
- A program which initiates the distribution of tracers
- and advects those tracers in a time evolving velocity field.
- Called and used from the CitCOM finite element code.
- Written 2/96 M. Gurnis for Citcom in cartesian geometry
- Modified by Lijie in 1998 and by Vlad and Eh in 2005 for CitcomS
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <mpi.h>
-#include <math.h>
-#include <sys/types.h>
-#ifdef HAVE_MALLOC_H
-#include <malloc.h>
-#endif
-#include <stdlib.h> /* for "system" command */
-
-#include "element_definitions.h"
-#include "global_defs.h"
-#include "parsing.h"
-
-void tracer_input(struct All_variables *E)
-{
- int m=E->parallel.me;
-
- input_int("tracer",&(E->control.tracer),"0",m);
- input_string("tracer_file",E->control.tracer_file,"tracer.dat",m);
-}
-
-
-void tracer_initial_settings(E)
- struct All_variables *E;
-{
- void tracer_setup();
- void tracer_advection();
-
- E->problem_tracer_setup=tracer_setup;
- E->problem_tracer_advection=tracer_advection;
-}
-
-
-void tracer_setup(E)
- struct All_variables *E;
-{
- int i,j,k,node;
- int m,ntr;
- float idummy,xdummy,ydummy,zdummy;
- FILE *fp;
- int nox,noy,noz,gnox,gnoy,gnoz;
- char output_file[255];
- MPI_Comm world;
- MPI_Status status;
-
- gnox=E->mesh.nox;
- gnoy=E->mesh.noy;
- gnoz=E->mesh.noz;
- nox=E->lmesh.nox;
- noy=E->lmesh.noy;
- noz=E->lmesh.noz;
-
- sprintf(output_file,"%s",E->control.tracer_file);
- fp=fopen(output_file,"r");
- if (fp == NULL) {
- fprintf(E->fp,"(Tracer_advection #1) Cannot open %s\n", output_file);
- exit(8);
- }
- fscanf(fp,"%d",&(E->Tracer.NUM_TRACERS));
-
- E->Tracer.tracer_x=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.tracer_y=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.tracer_z=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- E->Tracer.itcolor=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
-
- E->Tracer.x_space=(float*) malloc((gnox+1)*sizeof(float));
- E->Tracer.y_space=(float*) malloc((gnoy+1)*sizeof(float));
- E->Tracer.z_space=(float*) malloc((gnoz+1)*sizeof(float));
-
-
-
- /***comment by Vlad 03/15/2005
- each processor holds its own number of tracers
- ***/
-
- ntr=1;
- for(i=1;i<=E->Tracer.NUM_TRACERS;i++) {
- fscanf(fp,"%f %f %f %f", &idummy, &xdummy, &ydummy, &zdummy);
- if(xdummy >= E->sx[1][1][1] && xdummy <= E->sx[1][1][nox*noy*noz]) {
- if(ydummy >= E->sx[1][2][1] && ydummy <= E->sx[1][2][nox*noy*noz]) {
- if(zdummy >= E->sx[1][3][1] && zdummy <= E->sx[1][3][nox*noy*noz]) {
- E->Tracer.itcolor[ntr]=idummy;
- E->Tracer.tracer_x[ntr]=xdummy;
- E->Tracer.tracer_y[ntr]=ydummy;
- E->Tracer.tracer_z[ntr]=zdummy;
- ntr++;
- }
- }
- }
- }
-
- /***comment by Vlad 3/30/2005
- E->Tracer.LOCAL_NUM_TRACERS is the initial number
- of tracers in each processor
- ***/
-
- E->Tracer.LOCAL_NUM_TRACERS=ntr-1;
-
-
-
- /***comment by Vlad 1/26/2005
- reading the local mesh coordinate
- ***/
-
-
-
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(i=1;i<=nox;i++)
- {
- j=1;
- k=1;
- node=k+(i-1)*noz+(j-1)*nox*noz;
- E->Tracer.x_space[i]=E->sx[m][1][node];
- }
-
- for(j=1;j<=noy;j++)
- {
- i=1;
- k=1;
- node=k+(i-1)*noz+(j-1)*nox*noz;
- E->Tracer.y_space[j]=E->sx[m][2][node];
- }
-
- for(k=1;k<=noz;k++)
- {
- i=1;
- j=1;
- node=k+(i-1)*noz+(j-1)*nox*noz;
- E->Tracer.z_space[k]=E->sx[m][3][node];
- }
-
- } /* end of m */
-
-
- return;
-
-}
-
-void tracer_advection(E)
- struct All_variables *E;
-{
- int i,j,k,l,m,n,o,p;
- int n_x,n_y,n_z;
- int node1,node2,node3,node4,node5,node6,node7,node8;
- int nno,nox,noy,noz,gnox,gnoy,gnoz;
- int iteration;
- float x_tmp, y_tmp, z_tmp;
- float volume, tr_dx, tr_dy, tr_dz, dx, dy, dz;
- float w1,w2,w3,w4,w5,w6,w7,w8;
- float tr_v[NCS][4];
- MPI_Comm world;
- MPI_Status status[4];
- MPI_Status status_count;
- MPI_Request request[4];
- float xmin,xmax,ymin,ymax,zmin,zmax;
-
-
- float *tr_color_1,*tr_x_1,*tr_y_1,*tr_z_1;
- float *tr_color_new[13],*tr_x_new[13],*tr_y_new[13],*tr_z_new[13];
- int *Left_proc_list,*Right_proc_list;
- int *jump_new,*count_new;
- int *jj;
-
-
- int proc;
- int Previous_num_tracers,Current_num_tracers;
- int locx,locy,locz;
- int left,right,up,down,back,front;
-
-
-
- world=E->parallel.world;
-
-
- gnox=E->mesh.nox;
- gnoy=E->mesh.noy;
- gnoz=E->mesh.noz;
- nox=E->lmesh.nox;
- noy=E->lmesh.noy;
- noz=E->lmesh.noz;
- nno=nox*noy*noz;
-
- Left_proc_list=(int*) malloc(6*sizeof(int));
- Right_proc_list=(int*) malloc(6*sizeof(int));
- jump_new=(int*) malloc(6*sizeof(int));
- count_new=(int*) malloc(6*sizeof(int));
- jj=(int*) malloc(6*sizeof(int));
-
- tr_x_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_y_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_z_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_color_1=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
-
- for(i=0;i<=11;i++){
- tr_color_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_x_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_y_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- tr_z_new[i]=(float*) malloc((E->Tracer.NUM_TRACERS+1)*sizeof(float));
- }
-
-
- /*** comment by Vlad 3/25/2005
- This part of code gets the bounding box of the local mesh.
- ***/
-
- xmin=E->Tracer.x_space[1];
- xmax=E->Tracer.x_space[nox];
- ymin=E->Tracer.y_space[1];
- ymax=E->Tracer.y_space[noy];
- zmin=E->Tracer.z_space[1];
- zmax=E->Tracer.z_space[noz];
-
- /*fprintf(stderr,"%d %d\n", E->parallel.nprocx, E->parallel.loc2proc_map[0][0][1][0]);*/
-
- /*** comment by Tan2 1/25/2005
- Copy the velocity array.
- ***/
-
-
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(i=1;i<=nno;i++)
- for(j=1;j<=3;j++) {
- E->GV[m][j][i]=E->sphere.cap[m].V[j][i];
- }
- }
-
- /*** comment by vlad 03/17/2005
- advecting tracers in each processor
- ***/
-
-
- for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++) {
-
- n_x=0;
- n_y=0;
- n_z=0;
-
- /*printf("%s %d %s %d %s %f\n","La pasul",E->monitor.solution_cycles, "Procesorul", E->parallel.me, "advecteaza tracerul",E->Tracer.tracer_y[n]);
-
- //printf("%d %d %d %f %f %f %f\n", E->monitor.solution_cycles, E->parallel.me,E->Tracer.LOCAL_NUM_TRACERS,E->Tracer.itcolor[n],E->Tracer.tracer_x[n],E->Tracer.tracer_y[n],E->Tracer.tracer_z[n]);
- //printf("%d %d %d %f %f %f %f\n", E->monitor.solution_cycles, E->parallel.me,E->Tracer.LOCAL_NUM_TRACERS,tr_color[n],tr_x[n],tr_y[n],tr_z[n]);
- */
-
- /* mid point method uses 2 iterations */
-
- for(iteration=1;iteration<=2;iteration++)
- {
-
- if(iteration ==1) {
- x_tmp=E->Tracer.tracer_x[n];
- y_tmp=E->Tracer.tracer_y[n];
- z_tmp=E->Tracer.tracer_z[n];
- }
-
- /*** comment by Tan2 1/25/2005
- Find the element that contains the tracer.
-
- nodex n_x n_x+1
- | * |
- <----------->
- tr_dx
-
- <-------------------------------->
- dx
- ***/
-
- for(i=1;i<nox;i++) {
- if(x_tmp >= E->Tracer.x_space[i] && x_tmp <= E->Tracer.x_space[i+1]) {
- tr_dx=x_tmp-E->Tracer.x_space[i];
- dx=E->Tracer.x_space[i+1]-E->Tracer.x_space[i];
- n_x=i;
- }
- }
-
- for(j=1;j<noy;j++) {
- if(y_tmp >= E->Tracer.y_space[j] && y_tmp <= E->Tracer.y_space[j+1]) {
- tr_dy=y_tmp-E->Tracer.y_space[j];
- dy=E->Tracer.y_space[j+1]-E->Tracer.y_space[j];
- n_y=j;
- }
- }
-
- for(k=1;k<noz;k++) {
- if(z_tmp >= E->Tracer.z_space[k] && z_tmp <= E->Tracer.z_space[k+1]) {
- tr_dz=z_tmp-E->Tracer.z_space[k];
- dz=E->Tracer.z_space[k+1]-E->Tracer.z_space[k];
- n_z=k;
- }
- }
-
-
-
- /*** comment by Tan2 1/25/2005
- Calculate shape functions from tr_dx, tr_dy, tr_dz
- This assumes linear element
- ***/
-
-
- /* compute volumetic weighting functions */
- w1=tr_dx*tr_dz*tr_dy;
- w2=(dx-tr_dx)*tr_dz*tr_dy;
- w3=tr_dx*(dz-tr_dz)*tr_dy;
- w4=(dx-tr_dx)*(dz-tr_dz)*tr_dy;
- w5=tr_dx*tr_dz*(dy-tr_dy);
- w6=(dx-tr_dx)*tr_dz*(dy-tr_dy);
- w7=tr_dx*(dz-tr_dz)*(dy-tr_dy);
- w8=(dx-tr_dx)*(dz-tr_dz)*(dy-tr_dy);
-
-
- volume=dx*dz*dy;
-
-
- /*** comment by Tan2 1/25/2005
- Calculate the 8 node numbers of current element
- ***/
-
- node1 = n_z + (n_x-1)*noz + (n_y-1)*noz*nox;
- node2 = n_z + n_x*noz + (n_y-1)*noz*nox;
- node3 = n_z+1 + (n_x-1)*noz + (n_y-1)*noz*nox;
- node4 = n_z+1 + n_x*noz + (n_y-1)*noz*nox;
- node5 = n_z + (n_x-1)*noz + n_y*noz*nox;
- node6 = n_z + n_x*noz + n_y*noz*nox;
- node7 = n_z+1 + (n_x-1)*noz + n_y*noz*nox;
- node8 = n_z+1 + n_x*noz + n_y*noz*nox;
-
-
- /* printf("%d %d %d %d %d %d %d %d %d\n",E->parallel.me, node1,node2,node3,node4,node5,node6,node7,node8);
- //printf("%d %f %f %f %f %f %f %f %f\n", E->parallel.me, E->GV[1][2][node1], E->GV[1][2][node2], E->GV[1][2][node3], E->GV[1][2][node4], E->GV[1][2][node5], E->GV[1][2][node6], E->GV[1][2][node7], E->GV[1][2][node8]);
- */
-
- /*** comment by Tan2 1/25/2005
- Interpolate the velocity on the tracer position
- ***/
-
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(j=1;j<=3;j++) {
- tr_v[m][j]=w8*E->GV[m][j][node1]
- +w7*E->GV[m][j][node2]
- +w6*E->GV[m][j][node3]
- +w5*E->GV[m][j][node4]
- +w4*E->GV[m][j][node5]
- +w3*E->GV[m][j][node6]
- +w2*E->GV[m][j][node7]
- +w1*E->GV[m][j][node8];
- tr_v[m][j]=tr_v[m][j]/volume;
-
- }
-
- /*** comment by Tan2 1/25/2005
- advect tracer using mid-point method (2nd order accuracy)
- ***/
-
- /* mid point method */
-
- if(iteration == 1) {
- x_tmp = x_tmp + (E->advection.timestep/2.0)*tr_v[m][1]/E->Tracer.z_space[n_z];
- y_tmp = y_tmp + (E->advection.timestep/2.0)*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
- z_tmp = z_tmp + (E->advection.timestep/2.0)*tr_v[m][3];
- }
- if( iteration == 2) {
- E->Tracer.tracer_x[n] += E->advection.timestep*tr_v[m][1]/E->Tracer.z_space[n_z];
- E->Tracer.tracer_y[n] += E->advection.timestep*tr_v[m][2]/(E->Tracer.z_space[n_z]*sin(E->Tracer.x_space[n_x]));
- E->Tracer.tracer_z[n] += E->advection.timestep*tr_v[m][3];
-
- }
-
-
-
- } /* end of m */
-
-
- } /* end of iteration loop */
-
-
-
- }/* end of tracer loop */
-
-
- /*** Comment by Vlad 3/25/2005
- MPI for the tracer-advection code
- ***/
-
-
- m = 0;
-
- locx = E->parallel.me_loc[1];
- locy = E->parallel.me_loc[2];
- locz = E->parallel.me_loc[3];
-
- /* Am I the left-most proc.? If not, who is on my left? */
- if (locy == 0)
- left = -1;
- else
- left = E->parallel.loc2proc_map[m][locx][locy-1][locz];
-
- /* Am I the right-most proc.? If not, who is on my right? */
- if (locy == E->parallel.nprocy-1)
- right = -1;
- else
- right = E->parallel.loc2proc_map[m][locx][locy+1][locz];
-
- /* Am I the lower-most proc.? If not, who is beneath me? */
- if (locz == 0)
- down = -1;
- else
- down = E->parallel.loc2proc_map[m][locx][locy][locz-1];
-
- /* Am I the upper-most proc.? If not, who is above me? */
- if (locz == E->parallel.nprocz-1)
- up = -1;
- else
- up = E->parallel.loc2proc_map[m][locx][locy][locz+1];
-
- /* Am I the back-most proc.? If not, who is behind me? */
- if (locx == 0)
- back = -1;
- else
- back = E->parallel.loc2proc_map[m][locx-1][locy][locz];
-
- /* Am I the front-most proc.? If not, who is in front of me? */
- if (locx == E->parallel.nprocx-1)
- front = -1;
- else
- front = E->parallel.loc2proc_map[m][locx+1][locy][locz];
-
-
- Left_proc_list[0]=left;
- Left_proc_list[1]=right;
- Left_proc_list[2]=down;
- Left_proc_list[3]=up;
- Left_proc_list[4]=back;
- Left_proc_list[5]=front;
-
- Right_proc_list[0]=right;
- Right_proc_list[1]=left;
- Right_proc_list[2]=up;
- Right_proc_list[3]=down;
- Right_proc_list[4]=front;
- Right_proc_list[5]=back;
-
- jump_new[0]=0;
- jump_new[1]=0;
- jump_new[2]=0;
- jump_new[3]=0;
- jump_new[4]=0;
- jump_new[5]=0;
-
- count_new[0]=0;
- count_new[1]=0;
- count_new[2]=0;
- count_new[3]=0;
- count_new[4]=0;
- count_new[5]=0;
-
- jj[0]=1;
- jj[1]=0;
- jj[2]=3;
- jj[3]=2;
- jj[4]=5;
- jj[5]=4;
-
- E->Tracer.TEMP_TRACERS=0;
- Current_num_tracers=0;
-
- for(i=0;i<=11;i++){
- for(j=0;j<=E->Tracer.NUM_TRACERS;j++){
- tr_color_new[i][j]=999;
- tr_x_new[i][j]=999;
- tr_y_new[i][j]=999;
- tr_z_new[i][j]=999;
-
- tr_color_1[j]=999;
- tr_x_1[j]=999;
- tr_y_1[j]=999;
- tr_z_1[j]=999;
- }
- }
-
-
- i=0;
- j=0;
- k=0;
- l=0;
- m=0;
- o=0;
- p=0;
-
-
- for(n=1;n<=E->Tracer.LOCAL_NUM_TRACERS;n++){
-
- if(E->Tracer.tracer_y[n]>ymax) {
- if(E->Tracer.tracer_y[n]+100 != 100) {
- tr_color_new[0][i]=E->Tracer.itcolor[n];
- tr_x_new[0][i]=E->Tracer.tracer_x[n];
- tr_y_new[0][i]=E->Tracer.tracer_y[n];
- tr_z_new[0][i]=E->Tracer.tracer_z[n];
- i++;
- jump_new[0]=i;
- }
- }
- else if(E->Tracer.tracer_y[n]<ymin) {
- if(E->Tracer.tracer_y[n]+100 != 100) {
- tr_color_new[1][j]=E->Tracer.itcolor[n];
- tr_x_new[1][j]=E->Tracer.tracer_x[n];
- tr_y_new[1][j]=E->Tracer.tracer_y[n];
- tr_z_new[1][j]=E->Tracer.tracer_z[n];
- j++;
- jump_new[1]=j;
- }
- }
- else if(E->Tracer.tracer_z[n]>zmax) {
- if(E->Tracer.tracer_z[n]+100 != 100) {
- tr_color_new[2][k]=E->Tracer.itcolor[n];
- tr_x_new[2][k]=E->Tracer.tracer_x[n];
- tr_y_new[2][k]=E->Tracer.tracer_y[n];
- tr_z_new[2][k]=E->Tracer.tracer_z[n];
- k++;
- jump_new[2]=k;
- }
- }
- else if(E->Tracer.tracer_z[n]<zmin) {
- if(E->Tracer.tracer_z[n]+100 != 100) {
- tr_color_new[3][l]=E->Tracer.itcolor[n];
- tr_x_new[3][l]=E->Tracer.tracer_x[n];
- tr_y_new[3][l]=E->Tracer.tracer_y[n];
- tr_z_new[3][l]=E->Tracer.tracer_z[n];
- l++;
- jump_new[3]=l;
- }
- }
-
- else if(E->Tracer.tracer_x[n]>xmax) {
- if(E->Tracer.tracer_x[n]+100 != 100) {
- tr_color_new[4][m]=E->Tracer.itcolor[n];
- tr_x_new[4][m]=E->Tracer.tracer_x[n];
- tr_y_new[4][m]=E->Tracer.tracer_y[n];
- tr_z_new[4][m]=E->Tracer.tracer_z[n];
- m++;
- jump_new[4]=m;
- }
- }
- else if(E->Tracer.tracer_x[n]<xmin) {
- if(E->Tracer.tracer_x[n]+100 != 100) {
- tr_color_new[5][o]=E->Tracer.itcolor[n];
- tr_x_new[5][o]=E->Tracer.tracer_x[n];
- tr_y_new[5][o]=E->Tracer.tracer_y[n];
- tr_z_new[5][o]=E->Tracer.tracer_z[n];
- o++;
- jump_new[5]=o;
- }
- }
-
- else {
- tr_color_1[p]=E->Tracer.itcolor[n];
- tr_x_1[p]=E->Tracer.tracer_x[n];
- tr_y_1[p]=E->Tracer.tracer_y[n];
- tr_z_1[p]=E->Tracer.tracer_z[n];
- p++;
- }
- }
-
- Previous_num_tracers=E->Tracer.LOCAL_NUM_TRACERS;
- Current_num_tracers=Previous_num_tracers-jump_new[0]-jump_new[1]-jump_new[2]-jump_new[3]-jump_new[4]-jump_new[5];
-
- for(p=1;p<=Current_num_tracers;p++){
- E->Tracer.itcolor[p]=tr_color_1[p-1];
- E->Tracer.tracer_x[p]=tr_x_1[p-1];
- E->Tracer.tracer_y[p]=tr_y_1[p-1];
- E->Tracer.tracer_z[p]=tr_z_1[p-1];
- }
-
-
- for(i=0;i<=5;i++){
-
- j=jj[i];
-
- if (Left_proc_list[i] >= 0) {
- proc=Left_proc_list[i];
- MPI_Irecv(tr_color_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 11+i, world, &request[0]);
- MPI_Irecv(tr_x_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 12+i, world, &request[1]);
- MPI_Irecv(tr_y_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 13+i, world, &request[2]);
- MPI_Irecv(tr_z_new[i+6], E->Tracer.NUM_TRACERS, MPI_FLOAT, proc, 14+i, world, &request[3]);
- }
-
- if (Right_proc_list[i] >= 0) {
- proc=Right_proc_list[i];
- MPI_Send(tr_color_new[i], jump_new[i], MPI_FLOAT, proc, 11+i, world);
- MPI_Send(tr_x_new[i], jump_new[i], MPI_FLOAT, proc, 12+i, world);
- MPI_Send(tr_y_new[i], jump_new[i], MPI_FLOAT, proc, 13+i, world);
- MPI_Send(tr_z_new[i], jump_new[i], MPI_FLOAT, proc, 14+i, world);
- }
-
- if (Left_proc_list[i] >= 0) {
- MPI_Waitall(4, request, status);
- status_count = status[0];
- MPI_Get_count(&status_count, MPI_FLOAT, &count_new[i]);
- }
-
-
- E->Tracer.TEMP_TRACERS=E->Tracer.TEMP_TRACERS+count_new[i]-jump_new[i];
- E->Tracer.LOCAL_NUM_TRACERS=Previous_num_tracers+E->Tracer.TEMP_TRACERS;
-
-
-
- if(i <= 1){
- for(n=Current_num_tracers+count_new[j];n<=Current_num_tracers+count_new[i]+count_new[j]-1;n++) {
- m=Current_num_tracers+count_new[j];
- E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
- E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
- E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
- E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
-
- }
- }
-
-
- else if (i <= 4) {
- for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];n<=Current_num_tracers+count_new[0]+count_new[1]+count_new[i]+count_new[j]-1;n++) {
- m=Current_num_tracers+count_new[0]+count_new[1]+count_new[j];
- E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
- E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
- E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
- E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
-
- }
- }
-
- else {
- for(n=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];n<=E->Tracer.LOCAL_NUM_TRACERS-1;n++) {
- m=Current_num_tracers+count_new[0]+count_new[1]+count_new[2]+count_new[3]+count_new[j];
- E->Tracer.itcolor[n+1]=tr_color_new[i+6][n-m];
- E->Tracer.tracer_x[n+1]=tr_x_new[i+6][n-m];
- E->Tracer.tracer_y[n+1]=tr_y_new[i+6][n-m];
- E->Tracer.tracer_z[n+1]=tr_z_new[i+6][n-m];
-
- }
- }
-
-
- }
-
-
- free (tr_color_1);
- free (tr_x_1);
- free (tr_y_1);
- free (tr_z_1);
- for(i=0;i<=11;i++){
- free (tr_color_new[i]);
- free (tr_x_new[i]);
- free (tr_y_new[i]);
- free (tr_z_new[i]);
- }
-
-
- return;
-}
Added: mc/3D/CitcomS/branches/compressible/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Tracer_setup.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Tracer_setup.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -0,0 +1,67 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *
+ *<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>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+/*
+
+ Tracer_setup.c
+
+ A program which initiates the distribution of tracers
+ and advects those tracers in a time evolving velocity field.
+ Called and used from the CitCOM finite element code.
+ Written 2/96 M. Gurnis for Citcom in cartesian geometry
+ Modified by Lijie in 1998 and by Vlad and Eh in 2005 for CitcomS
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "global_defs.h"
+#include "parsing.h"
+
+
+void tracer_input(struct All_variables *E)
+{
+ int m=E->parallel.me;
+
+ input_int("tracer",&(E->control.tracer),"0",m);
+ input_string("tracer_file",E->control.tracer_file,"tracer.dat",m);
+}
+
+
+void tracer_initial_settings(E)
+ struct All_variables *E;
+{
+ void regional_tracer_setup();
+ void regional_tracer_advection();
+
+ if(E->parallel.nprocxy == 1) {
+ E->problem_tracer_setup = regional_tracer_setup;
+ E->problem_tracer_advection = regional_tracer_advection;
+ }
+}
Modified: mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -36,69 +36,81 @@
#include "global_defs.h"
#include "parsing.h"
+static void apply_low_visc_wedge_channel(struct All_variables *E, float **evisc);
+static void low_viscosity_channel_factor(struct All_variables *E, float *F);
+static void low_viscosity_wedge_factor(struct All_variables *E, float *F);
+
void viscosity_system_input(struct All_variables *E)
{
- int m=E->parallel.me;
- int i;
+ int m=E->parallel.me;
+ int i;
- /* default values .... */
- for(i=0;i<40;i++) {
- E->viscosity.N0[i]=1.0;
- E->viscosity.T[i] = 0.0;
- E->viscosity.Z[i] = 0.0;
- E->viscosity.E[i] = 0.0;
- }
+ /* default values .... */
+ for(i=0;i<40;i++) {
+ E->viscosity.N0[i]=1.0;
+ E->viscosity.T[i] = 0.0;
+ E->viscosity.Z[i] = 0.0;
+ E->viscosity.E[i] = 0.0;
+ }
- /* read in information */
- input_boolean("VISC_UPDATE",&(E->viscosity.update_allowed),"on",m);
- input_int("rheol",&(E->viscosity.RHEOL),"3",m);
- input_int("num_mat",&(E->viscosity.num_mat),"1",m);
- input_float_vector("visc0",E->viscosity.num_mat,(E->viscosity.N0),m);
+ /* read in information */
+ input_boolean("VISC_UPDATE",&(E->viscosity.update_allowed),"on",m);
+ input_int("rheol",&(E->viscosity.RHEOL),"3",m);
+ input_int("num_mat",&(E->viscosity.num_mat),"1",m);
+ input_float_vector("visc0",E->viscosity.num_mat,(E->viscosity.N0),m);
- input_boolean("TDEPV",&(E->viscosity.TDEPV),"on",m);
- if (E->viscosity.TDEPV) {
- input_float_vector("viscT",E->viscosity.num_mat,(E->viscosity.T),m);
- input_float_vector("viscE",E->viscosity.num_mat,(E->viscosity.E),m);
- input_float_vector("viscZ",E->viscosity.num_mat,(E->viscosity.Z),m);
- }
+ input_boolean("TDEPV",&(E->viscosity.TDEPV),"on",m);
+ if (E->viscosity.TDEPV) {
+ input_float_vector("viscT",E->viscosity.num_mat,(E->viscosity.T),m);
+ input_float_vector("viscE",E->viscosity.num_mat,(E->viscosity.E),m);
+ input_float_vector("viscZ",E->viscosity.num_mat,(E->viscosity.Z),m);
+ }
- E->viscosity.sdepv_misfit = 1.0;
- input_boolean("SDEPV",&(E->viscosity.SDEPV),"off",m);
- if (E->viscosity.SDEPV) {
- input_float("sdepv_misfit",&(E->viscosity.sdepv_misfit),"0.001",m);
- input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
- }
+ E->viscosity.sdepv_misfit = 1.0;
+ input_boolean("SDEPV",&(E->viscosity.SDEPV),"off",m);
+ if (E->viscosity.SDEPV) {
+ input_float("sdepv_misfit",&(E->viscosity.sdepv_misfit),"0.001",m);
+ input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
+ }
- input_boolean("VMAX",&(E->viscosity.MAX),"off",m);
- if (E->viscosity.MAX)
- input_float("visc_max",&(E->viscosity.max_value),"1e22,1,nomax",m);
+ input_boolean("low_visc_channel",&(E->viscosity.channel),"off",m);
+ input_boolean("low_visc_wedge",&(E->viscosity.wedge),"off",m);
- input_boolean("VMIN",&(E->viscosity.MIN),"off",m);
- if (E->viscosity.MIN)
- input_float("visc_min",&(E->viscosity.min_value),"1e20",m);
+ input_float("lv_min_radius",&(E->viscosity.lv_min_radius),"0.9764",m);
+ input_float("lv_max_radius",&(E->viscosity.lv_max_radius),"0.9921",m);
+ input_float("lv_channel_thickness",&(E->viscosity.lv_channel_thickness),"0.0047",m);
+ input_float("lv_reduction",&(E->viscosity.lv_reduction),"0.5",m);
- return;
+ input_boolean("VMAX",&(E->viscosity.MAX),"off",m);
+ if (E->viscosity.MAX)
+ input_float("visc_max",&(E->viscosity.max_value),"1e22,1,nomax",m);
+
+ input_boolean("VMIN",&(E->viscosity.MIN),"off",m);
+ if (E->viscosity.MIN)
+ input_float("visc_min",&(E->viscosity.min_value),"1e20",m);
+
+ return;
}
void viscosity_input(struct All_variables *E)
{
- int m = E->parallel.me;
+ int m = E->parallel.me;
- input_string("Viscosity",E->viscosity.STRUCTURE,NULL,m);
- input_int ("visc_smooth_method",&(E->viscosity.smooth_cycles),"0",m);
+ input_string("Viscosity",E->viscosity.STRUCTURE,NULL,m);
+ input_int ("visc_smooth_method",&(E->viscosity.smooth_cycles),"0",m);
- if ( strcmp(E->viscosity.STRUCTURE,"system") == 0)
- E->viscosity.FROM_SYSTEM = 1;
- else
- E->viscosity.FROM_SYSTEM = 0;
+ if ( strcmp(E->viscosity.STRUCTURE,"system") == 0)
+ E->viscosity.FROM_SYSTEM = 1;
+ else
+ E->viscosity.FROM_SYSTEM = 0;
- if (E->viscosity.FROM_SYSTEM)
- viscosity_system_input(E);
+ if (E->viscosity.FROM_SYSTEM)
+ viscosity_system_input(E);
- return;
+ return;
}
@@ -126,66 +138,72 @@
const int vpts = vpoints[E->mesh.nsd];
if(E->viscosity.TDEPV)
- visc_from_T(E,evisc,propogate);
+ visc_from_T(E,evisc,propogate);
else
- visc_from_mat(E,evisc);
+ visc_from_mat(E,evisc);
if(E->viscosity.SDEPV)
- visc_from_S(E,evisc,propogate);
+ visc_from_S(E,evisc,propogate);
+
+ apply_low_visc_wedge_channel(E, evisc);
+
+
+ /* min/max cut-off */
+
if(E->viscosity.MAX) {
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nel;i++)
- for(j=1;j<=vpts;j++)
- if(evisc[m][(i-1)*vpts + j] > E->viscosity.max_value)
- evisc[m][(i-1)*vpts + j] = E->viscosity.max_value;
- }
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.nel;i++)
+ for(j=1;j<=vpts;j++)
+ if(evisc[m][(i-1)*vpts + j] > E->viscosity.max_value)
+ evisc[m][(i-1)*vpts + j] = E->viscosity.max_value;
+ }
if(E->viscosity.MIN) {
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nel;i++)
- for(j=1;j<=vpts;j++)
- if(evisc[m][(i-1)*vpts + j] < E->viscosity.min_value)
- evisc[m][(i-1)*vpts + j] = E->viscosity.min_value;
- }
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.nel;i++)
+ for(j=1;j<=vpts;j++)
+ if(evisc[m][(i-1)*vpts + j] < E->viscosity.min_value)
+ evisc[m][(i-1)*vpts + j] = E->viscosity.min_value;
+ }
/*
- if (E->control.verbose) {
- fprintf(E->fp_out,"output_evisc \n");
- for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ if (E->control.verbose) {
+ fprintf(E->fp_out,"output_evisc \n");
+ for(m=1;m<=E->sphere.caps_per_proc;m++) {
fprintf(E->fp_out,"output_evisc for cap %d\n",E->sphere.capid[m]);
for(i=1;i<=E->lmesh.nel;i++)
- fprintf(E->fp_out,"%d %d %f %f\n",i,E->mat[m][i],evisc[m][(i-1)*vpts+1],evisc[m][(i-1)*vpts+7]);
+ fprintf(E->fp_out,"%d %d %f %f\n",i,E->mat[m][i],evisc[m][(i-1)*vpts+1],evisc[m][(i-1)*vpts+7]);
}
- fflush(E->fp_out);
- }
+ fflush(E->fp_out);
+ }
*/
visc_from_gint_to_nodes(E,evisc,visc,E->mesh.levmax);
visc_from_nodes_to_gint(E,visc,evisc,E->mesh.levmax);
-/* visc_to_node_interpolate(E,evisc,visc);
-*/
+ /* visc_to_node_interpolate(E,evisc,visc);
+ */
-/* for(m=1;m<=E->sphere.caps_per_proc;m++) {
- for(i=1;i<=E->lmesh.nel;i++)
- if (i%E->lmesh.elz==0) {
+ /* for(m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(i=1;i<=E->lmesh.nel;i++)
+ if (i%E->lmesh.elz==0) {
fprintf(E->fp_out,"%.4e %.4e %.4e %5d %2d\n",E->eco[m][i].centre[1],E->eco[m][i].centre[2],log10(evisc[m][(i-1)*vpts+1]),i,E->mat[m][i]);
}
- } */
- return;
+ } */
+ return;
}
void initial_viscosity(struct All_variables *E)
{
- if (E->viscosity.FROM_SYSTEM)
- get_system_viscosity(E,1,E->EVI[E->mesh.levmax],E->VI[E->mesh.levmax]);
+ if (E->viscosity.FROM_SYSTEM)
+ get_system_viscosity(E,1,E->EVI[E->mesh.levmax],E->VI[E->mesh.levmax]);
- return;
+ return;
}
@@ -196,13 +214,13 @@
int i,m,jj;
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nel;i++)
- for(jj=1;jj<=vpoints[E->mesh.nsd];jj++)
- EEta[m][ (i-1)*vpoints[E->mesh.nsd]+jj ]=E->viscosity.N0[E->mat[m][i]-1];
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.nel;i++)
+ for(jj=1;jj<=vpoints[E->mesh.nsd];jj++)
+ EEta[m][ (i-1)*vpoints[E->mesh.nsd]+jj ]=E->viscosity.N0[E->mat[m][i]-1];
return;
- }
+}
void visc_from_T(E,EEta,propogate)
struct All_variables *E;
@@ -224,188 +242,188 @@
switch (E->viscosity.RHEOL) {
case 1:
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=nel;i++) {
- l = E->mat[m][i];
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=nel;i++) {
+ l = E->mat[m][i];
- if(E->control.mat_control==0)
- tempa = E->viscosity.N0[l-1];
- else if(E->control.mat_control==1)
- tempa = E->viscosity.N0[l-1]*E->VIP[m][i];
+ if(E->control.mat_control==0)
+ tempa = E->viscosity.N0[l-1];
+ else if(E->control.mat_control==1)
+ tempa = E->viscosity.N0[l-1]*E->VIP[m][i];
- for(kk=1;kk<=ends;kk++) {
- TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- }
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk] = E->T[m][E->ien[m][i].node[kk]];
+ }
- for(jj=1;jj<=vpts;jj++) {
- temp=0.0;
- for(kk=1;kk<=ends;kk++) {
- temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
- }
+ for(jj=1;jj<=vpts;jj++) {
+ temp=0.0;
+ for(kk=1;kk<=ends;kk++) {
+ temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ }
- EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( E->viscosity.E[l-1] * (E->viscosity.T[l-1] - temp));
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( E->viscosity.E[l-1] * (E->viscosity.T[l-1] - temp));
- }
- }
- break;
+ }
+ }
+ break;
case 2:
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=nel;i++) {
- l = E->mat[m][i];
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=nel;i++) {
+ l = E->mat[m][i];
- if(E->control.mat_control==0)
- tempa = E->viscosity.N0[l-1];
- else if(E->control.mat_control==1)
- tempa = E->viscosity.N0[l-1]*E->VIP[m][i];
+ if(E->control.mat_control==0)
+ tempa = E->viscosity.N0[l-1];
+ else if(E->control.mat_control==1)
+ tempa = E->viscosity.N0[l-1]*E->VIP[m][i];
- for(kk=1;kk<=ends;kk++) {
- TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- }
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk] = E->T[m][E->ien[m][i].node[kk]];
+ }
- for(jj=1;jj<=vpts;jj++) {
- temp=0.0;
- for(kk=1;kk<=ends;kk++) {
- temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
- }
+ for(jj=1;jj<=vpts;jj++) {
+ temp=0.0;
+ for(kk=1;kk<=ends;kk++) {
+ temp += TT[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ }
- EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( -temp / E->viscosity.T[l-1]);
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( -temp / E->viscosity.T[l-1]);
- }
- }
- break;
+ }
+ }
+ break;
case 3:
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=nel;i++) {
- l = E->mat[m][i];
- tempa = E->viscosity.N0[l-1];
- j = 0;
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=nel;i++) {
+ l = E->mat[m][i];
+ tempa = E->viscosity.N0[l-1];
+ j = 0;
- for(kk=1;kk<=ends;kk++) {
- TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
- }
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk] = E->T[m][E->ien[m][i].node[kk]];
+ zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
+ }
- for(jj=1;jj<=vpts;jj++) {
- temp=0.0;
- zzz=0.0;
- for(kk=1;kk<=ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
- zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
- }
+ for(jj=1;jj<=vpts;jj++) {
+ temp=0.0;
+ zzz=0.0;
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk]=max(TT[kk],zero);
+ temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ }
- if(E->control.mat_control==0)
- EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+ if(E->control.mat_control==0)
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+ - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
- if(E->control.mat_control==1)
- EEta[m][ (i-1)*vpts + jj ] = tempa*E->VIP[m][i]*
- exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
- }
- }
- break;
+ if(E->control.mat_control==1)
+ EEta[m][ (i-1)*vpts + jj ] = tempa*E->VIP[m][i]*
+ exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+ - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+ }
+ }
+ break;
case 4:
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=nel;i++) {
- l = E->mat[m][i];
- tempa = E->viscosity.N0[l-1];
- j = 0;
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=nel;i++) {
+ l = E->mat[m][i];
+ tempa = E->viscosity.N0[l-1];
+ j = 0;
- for(kk=1;kk<=ends;kk++) {
- TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
- }
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk] = E->T[m][E->ien[m][i].node[kk]];
+ zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
+ }
- for(jj=1;jj<=vpts;jj++) {
- temp=0.0;
- zzz=0.0;
- for(kk=1;kk<=ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
- zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
- }
+ for(jj=1;jj<=vpts;jj++) {
+ temp=0.0;
+ zzz=0.0;
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk]=max(TT[kk],zero);
+ temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ }
-/* The viscosity formulation (dimensional) is: visc=visc0*exp[(Ea+p*Va)/R*T]
- Typical values for dry upper mantle are: Ea = 300 KJ/mol ; Va = 1.e-5 m^3/mol
- T=T0+DT*T'; where DT - temperature contrast (from Rayleigh number)
- T' - nondimensional temperature; T0 - surface tempereture (273 K)
- T=DT*[(T0/DT) + T'] => visc=visc0*exp{(Ea+p*Va)/R*DT*[(T0/DT) + T']}
- visc=visc0*exp{[(Ea/R*DT) + (p*Va/R*DT)]/[(T0/DT) + T']}
- so: E->viscosity.E = Ea/R*DT ; E->viscosity.Z = Va/R*DT
- p = zzz and E->viscosity.T = T0/DT */
+ /* The viscosity formulation (dimensional) is: visc=visc0*exp[(Ea+p*Va)/R*T]
+ Typical values for dry upper mantle are: Ea = 300 KJ/mol ; Va = 1.e-5 m^3/mol
+ T=T0+DT*T'; where DT - temperature contrast (from Rayleigh number)
+ T' - nondimensional temperature; T0 - surface tempereture (273 K)
+ T=DT*[(T0/DT) + T'] => visc=visc0*exp{(Ea+p*Va)/R*DT*[(T0/DT) + T']}
+ visc=visc0*exp{[(Ea/R*DT) + (p*Va/R*DT)]/[(T0/DT) + T']}
+ so: E->viscosity.E = Ea/R*DT ; E->viscosity.Z = Va/R*DT
+ p = zzz and E->viscosity.T = T0/DT */
- if(E->control.mat_control==0)
- EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
- / (E->viscosity.T[l-1]+temp) );
+ if(E->control.mat_control==0)
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
+ / (E->viscosity.T[l-1]+temp) );
- if(E->control.mat_control==1)
- EEta[m][ (i-1)*vpts + jj ] = tempa*E->VIP[m][i]*
- exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
- / (E->viscosity.T[l-1]+temp) );
+ if(E->control.mat_control==1)
+ EEta[m][ (i-1)*vpts + jj ] = tempa*E->VIP[m][i]*
+ exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
+ / (E->viscosity.T[l-1]+temp) );
- }
- }
- break;
+ }
+ }
+ break;
case 5:
- /* same as rheol 3, except alternative margin, VIP, formulation */
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=nel;i++) {
- l = E->mat[m][i];
- tempa = E->viscosity.N0[l-1];
- j = 0;
+ /* same as rheol 3, except alternative margin, VIP, formulation */
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=nel;i++) {
+ l = E->mat[m][i];
+ tempa = E->viscosity.N0[l-1];
+ j = 0;
- for(kk=1;kk<=ends;kk++) {
- TT[kk] = E->T[m][E->ien[m][i].node[kk]];
- zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
- }
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk] = E->T[m][E->ien[m][i].node[kk]];
+ zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
+ }
- for(jj=1;jj<=vpts;jj++) {
- temp=0.0;
- zzz=0.0;
- for(kk=1;kk<=ends;kk++) {
- TT[kk]=max(TT[kk],zero);
- temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
- zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
- }
+ for(jj=1;jj<=vpts;jj++) {
+ temp=0.0;
+ zzz=0.0;
+ for(kk=1;kk<=ends;kk++) {
+ TT[kk]=max(TT[kk],zero);
+ temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+ zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+ }
- if(E->control.mat_control==0)
- EEta[m][ (i-1)*vpts + jj ] = tempa*
- exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+ if(E->control.mat_control==0)
+ EEta[m][ (i-1)*vpts + jj ] = tempa*
+ exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+ - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
- if(E->control.mat_control==1) {
- visc1 = E->VIP[m][i];
- visc2 = 2.0/(1./visc1 + 1.);
- tempa_exp = tempa*
- exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
- - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
- visc1 = tempa*E->viscosity.max_value;
- if(tempa_exp > visc1) tempa_exp=visc1;
- EEta[m][ (i-1)*vpts + jj ] = visc2*tempa_exp;
- /* if(E->parallel.me == 0 && visc1 < 1.0e-03)
- fprintf(stderr,"%f %f %e %e %e\n",zzz,temp,visc1,visc2,
- EEta[m][ (i-1)*vpts + jj ]); */
- }
+ if(E->control.mat_control==1) {
+ visc1 = E->VIP[m][i];
+ visc2 = 2.0/(1./visc1 + 1.);
+ tempa_exp = tempa*
+ exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+ - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+ visc1 = tempa*E->viscosity.max_value;
+ if(tempa_exp > visc1) tempa_exp=visc1;
+ EEta[m][ (i-1)*vpts + jj ] = visc2*tempa_exp;
+ /* if(E->parallel.me == 0 && visc1 < 1.0e-03)
+ fprintf(stderr,"%f %f %e %e %e\n",zzz,temp,visc1,visc2,
+ EEta[m][ (i-1)*vpts + jj ]); */
+ }
- }
- }
- break;
+ }
+ }
+ break;
@@ -435,14 +453,14 @@
two = 2.0;
for(m=1;m<=E->sphere.caps_per_proc;m++) {
- strain_rate_2_inv(E,m,eedot,1);
+ strain_rate_2_inv(E,m,eedot,1);
- for(e=1;e<=nel;e++) {
- exponent1= one/E->viscosity.sdepv_expt[E->mat[m][e]-1];
- scale=pow(eedot[e],exponent1-one);
- for(jj=1;jj<=vpts;jj++)
- EEta[m][(e-1)*vpts + jj] = scale*pow(EEta[m][(e-1)*vpts+jj],exponent1);
- }
+ for(e=1;e<=nel;e++) {
+ exponent1= one/E->viscosity.sdepv_expt[E->mat[m][e]-1];
+ scale=pow(eedot[e],exponent1-one);
+ for(jj=1;jj<=vpts;jj++)
+ EEta[m][(e-1)*vpts + jj] = scale*pow(EEta[m][(e-1)*vpts+jj],exponent1);
+ }
}
free ((void *)eedot);
@@ -479,34 +497,34 @@
for(e=1;e<=nel;e++) {
- get_global_shape_fn(E,e,&GN,&GNx,&dOmega,2,sphere_key,rtf,lev,m);
+ get_global_shape_fn(E,e,&GN,&GNx,&dOmega,2,sphere_key,rtf,lev,m);
- velo_from_element(E,VV,m,e,sphere_key);
+ velo_from_element(E,VV,m,e,sphere_key);
- for(p=1;p<=dims;p++)
- for(q=1;q<=dims;q++)
- dudx[p][q] = 0.0;
+ for(p=1;p<=dims;p++)
+ for(q=1;q<=dims;q++)
+ dudx[p][q] = 0.0;
- for(i=1;i<=ends;i++)
+ for(i=1;i<=ends;i++)
+ for(p=1;p<=dims;p++)
+ for(q=1;q<=dims;q++)
+ dudx[p][q] += VV[p][i] * GNx.ppt[GNPXINDEX(q-1,i,1)];
+
for(p=1;p<=dims;p++)
- for(q=1;q<=dims;q++)
- dudx[p][q] += VV[p][i] * GNx.ppt[GNPXINDEX(q-1,i,1)];
+ for(q=1;q<=dims;q++)
+ edot[p][q] = dudx[p][q] + dudx[q][p];
- for(p=1;p<=dims;p++)
- for(q=1;q<=dims;q++)
- edot[p][q] = dudx[p][q] + dudx[q][p];
+ if (dims==2)
+ EEDOT[e] = edot[1][1]*edot[1][1] + edot[2][2]*edot[2][2]
+ + edot[1][2]*edot[1][2]*2.0;
- if (dims==2)
- EEDOT[e] = edot[1][1]*edot[1][1] + edot[2][2]*edot[2][2]
- + edot[1][2]*edot[1][2]*2.0;
+ else if (dims==3)
+ EEDOT[e] = edot[1][1]*edot[1][1] + edot[1][2]*edot[1][2]*2.0
+ + edot[2][2]*edot[2][2] + edot[2][3]*edot[2][3]*2.0
+ + edot[3][3]*edot[3][3] + edot[1][3]*edot[1][3]*2.0;
- else if (dims==3)
- EEDOT[e] = edot[1][1]*edot[1][1] + edot[1][2]*edot[1][2]*2.0
- + edot[2][2]*edot[2][2] + edot[2][3]*edot[2][3]*2.0
- + edot[3][3]*edot[3][3] + edot[1][3]*edot[1][3]*2.0;
+ }
- }
-
if(SQRT)
for(e=1;e<=nel;e++)
EEDOT[e] = sqrt(0.5 *EEDOT[e]);
@@ -520,70 +538,180 @@
void visc_to_node_interpolate(E,evisc,visc)
- struct All_variables *E;
- float **evisc,**visc;
+ struct All_variables *E;
+ float **evisc,**visc;
{
-/* void exchange_node_f(); */
-/* void get_global_shape_fn(); */
-/* void return_horiz_ave_f(); */
-/* void sphere_interpolate(); */
-/* void print_interpolated(); */
-/* void gather_TG_to_me0(); */
-/* void parallel_process_termination(); */
-/* int i,j,k,e,node,snode,m,nel2; */
-/* FILE *fp; */
-/* char output_file[255]; */
+ /* void exchange_node_f(); */
+ /* void get_global_shape_fn(); */
+ /* void return_horiz_ave_f(); */
+ /* void sphere_interpolate(); */
+ /* void print_interpolated(); */
+ /* void gather_TG_to_me0(); */
+ /* void parallel_process_termination(); */
+ /* int i,j,k,e,node,snode,m,nel2; */
+ /* FILE *fp; */
+ /* char output_file[255]; */
-/* float *TG,t,f,rad, Szz; */
+ /* float *TG,t,f,rad, Szz; */
-/* double time1,CPU_time0(),tww[9],rtf[4][9]; */
+ /* double time1,CPU_time0(),tww[9],rtf[4][9]; */
-/* struct Shape_function GN; */
-/* struct Shape_function_dA dOmega; */
-/* struct Shape_function_dx GNx; */
+ /* struct Shape_function GN; */
+ /* struct Shape_function_dA dOmega; */
+ /* struct Shape_function_dx GNx; */
-/* const int dims=E->mesh.nsd,dofs=E->mesh.dof; */
-/* const int vpts=vpoints[dims]; */
-/* const int ppts=ppoints[dims]; */
-/* const int ends=enodes[dims]; */
-/* const int nno=E->lmesh.nno; */
-/* const int lev=E->mesh.levmax; */
+ /* const int dims=E->mesh.nsd,dofs=E->mesh.dof; */
+ /* const int vpts=vpoints[dims]; */
+ /* const int ppts=ppoints[dims]; */
+ /* const int ends=enodes[dims]; */
+ /* const int nno=E->lmesh.nno; */
+ /* const int lev=E->mesh.levmax; */
-/* TG =(float *)malloc((E->sphere.nsf+1)*sizeof(float)); */
-/* for (i=E->sphere.nox;i>=1;i--) */
-/* for (j=1;j<=E->sphere.noy;j++) { */
-/* node = i + (j-1)*E->sphere.nox; */
-/* TG[node] = 0.0; */
-/* m = E->sphere.int_cap[node]; */
-/* e = E->sphere.int_ele[node]; */
+ /* TG =(float *)malloc((E->sphere.nsf+1)*sizeof(float)); */
+ /* for (i=E->sphere.nox;i>=1;i--) */
+ /* for (j=1;j<=E->sphere.noy;j++) { */
+ /* node = i + (j-1)*E->sphere.nox; */
+ /* TG[node] = 0.0; */
+ /* m = E->sphere.int_cap[node]; */
+ /* e = E->sphere.int_ele[node]; */
-/* if (m>0 && e>0) { */
-/* e=e+E->lmesh.elz-1; */
-/* TG[node] = log10(evisc[m][(e-1)*vpts+1]); */
-/* } */
-/* } */
+ /* if (m>0 && e>0) { */
+ /* e=e+E->lmesh.elz-1; */
+ /* TG[node] = log10(evisc[m][(e-1)*vpts+1]); */
+ /* } */
+ /* } */
-/* gather_TG_to_me0(E,TG); */
+ /* gather_TG_to_me0(E,TG); */
-/* if (E->parallel.me==E->parallel.nprocz-1) { */
-/* sprintf(output_file,"%s.evisc_intp",E->control.data_file); */
-/* fp=fopen(output_file,"w"); */
+ /* if (E->parallel.me==E->parallel.nprocz-1) { */
+ /* sprintf(output_file,"%s.evisc_intp",E->control.data_file); */
+ /* fp=fopen(output_file,"w"); */
-/* rad = 180/M_PI; */
-/* for (i=E->sphere.nox;i>=1;i--) */
-/* for (j=1;j<=E->sphere.noy;j++) { */
-/* node = i + (j-1)*E->sphere.nox; */
-/* t = 90-E->sphere.sx[1][node]*rad; */
-/* f = E->sphere.sx[2][node]*rad; */
-/* fprintf (fp,"%.3e %.3e %.4e\n",f,t,TG[node]); */
-/* } */
-/* fclose(fp); */
-/* } */
+ /* rad = 180/M_PI; */
+ /* for (i=E->sphere.nox;i>=1;i--) */
+ /* for (j=1;j<=E->sphere.noy;j++) { */
+ /* node = i + (j-1)*E->sphere.nox; */
+ /* t = 90-E->sphere.sx[1][node]*rad; */
+ /* f = E->sphere.sx[2][node]*rad; */
+ /* fprintf (fp,"%.3e %.3e %.4e\n",f,t,TG[node]); */
+ /* } */
+ /* fclose(fp); */
+ /* } */
-/* free((void *)TG); */
+ /* free((void *)TG); */
- return;
- }
+ return;
+}
+
+static void apply_low_visc_wedge_channel(struct All_variables *E, float **evisc)
+{
+ int i,j,jj,k,m,n,mm,ii,nn;
+ float temp1,temp2,*vvvis;
+
+ float THETA_LOCAL_ELEM_T, THETA_LOCAL_ELEM, FI_LOCAL_ELEM_T, FI_LOCAL_ELEM, R_LOCAL_ELEM_T, R_LOCAL_ELEM;
+ const int vpts = vpoints[E->mesh.nsd];
+
+ float *F;
+
+ F = (float *)malloc((E->lmesh.nel+1)*sizeof(float));
+ for(i=1 ; i<=E->lmesh.nel ; i++)
+ F[i] = 0.0;
+
+ /* if low viscosity channel ... */
+ if(E->viscosity.channel)
+ low_viscosity_channel_factor(E, F);
+
+
+ /* if low viscosity wedge ... */
+ if(E->viscosity.wedge)
+ low_viscosity_wedge_factor(E, F);
+
+
+ for(i=1 ; i<=E->lmesh.nel ; i++) {
+ if (F[i] != 0.0)
+ for(m = 1 ; m <= E->sphere.caps_per_proc ; m++) {
+ for(j=1;j<=vpts;j++) {
+ evisc[m][(i-1)*vpts + j] = F[i];
+ }
+ }
+ }
+
+
+ free(F);
+
+ return;
+}
+
+
+
+
+static void low_viscosity_channel_factor(struct All_variables *E, float *F)
+{
+ int i, n;
+ float THETA_LOCAL_ELEM_T, THETA_LOCAL_ELEM, FI_LOCAL_ELEM_T, FI_LOCAL_ELEM, R_LOCAL_ELEM_T, R_LOCAL_ELEM;
+
+
+ for(i=1 ; i<=E->lmesh.nel ; i++) {
+ THETA_LOCAL_ELEM = E->Tracer.THETA_LOC_ELEM[i];
+ FI_LOCAL_ELEM = E->Tracer.FI_LOC_ELEM[i];
+ R_LOCAL_ELEM = E->Tracer.R_LOC_ELEM[i];
+
+ if (R_LOCAL_ELEM >= E->viscosity.lv_min_radius && R_LOCAL_ELEM <= E->viscosity.lv_max_radius) {
+
+ for(n=1 ; n<=E->Tracer.LOCAL_NUM_TRACERS ; n++) {
+ if (E->Tracer.tracer_z[n] <= E->viscosity.lv_max_radius && E->Tracer.tracer_z[n] >= E->viscosity.lv_min_radius) {
+ THETA_LOCAL_ELEM_T = E->Tracer.THETA_LOC_ELEM_T[n];
+ FI_LOCAL_ELEM_T = E->Tracer.FI_LOC_ELEM_T[n];
+ R_LOCAL_ELEM_T = E->Tracer.R_LOC_ELEM_T[n];
+
+ if((R_LOCAL_ELEM >= R_LOCAL_ELEM_T) &&
+ (R_LOCAL_ELEM <= R_LOCAL_ELEM_T+E->viscosity.lv_channel_thickness) &&
+ (FI_LOCAL_ELEM == FI_LOCAL_ELEM_T) &&
+ (THETA_LOCAL_ELEM == THETA_LOCAL_ELEM_T)) {
+
+ F[i] = E->viscosity.lv_reduction;
+
+ }
+ }
+
+ }
+ }
+ }
+
+}
+
+
+
+static void low_viscosity_wedge_factor(struct All_variables *E, float *F)
+{
+ int i, n;
+ float THETA_LOCAL_ELEM_T, THETA_LOCAL_ELEM, FI_LOCAL_ELEM_T, FI_LOCAL_ELEM, R_LOCAL_ELEM_T, R_LOCAL_ELEM;
+
+ for(i=1 ; i<=E->lmesh.nel ; i++) {
+ THETA_LOCAL_ELEM = E->Tracer.THETA_LOC_ELEM[i];
+ FI_LOCAL_ELEM = E->Tracer.FI_LOC_ELEM[i];
+ R_LOCAL_ELEM = E->Tracer.R_LOC_ELEM[i];
+
+ if (R_LOCAL_ELEM <= E->viscosity.lv_max_radius && R_LOCAL_ELEM >= E->viscosity.lv_min_radius) {
+
+ for(n=1 ; n<=E->Tracer.LOCAL_NUM_TRACERS ; n++) {
+ if (E->Tracer.tracer_z[n] <= E->viscosity.lv_max_radius && E->Tracer.tracer_z[n] >= E->viscosity.lv_min_radius) {
+ THETA_LOCAL_ELEM_T = E->Tracer.THETA_LOC_ELEM_T[n];
+ FI_LOCAL_ELEM_T = E->Tracer.FI_LOC_ELEM_T[n];
+ R_LOCAL_ELEM_T = E->Tracer.R_LOC_ELEM_T[n];
+
+
+ if((R_LOCAL_ELEM >= R_LOCAL_ELEM_T) &&
+ (FI_LOCAL_ELEM == FI_LOCAL_ELEM_T) &&
+ (THETA_LOCAL_ELEM == THETA_LOCAL_ELEM_T)) {
+
+ F[i] = E->viscosity.lv_reduction;
+ }
+ }
+ }
+ }
+ }
+}
Modified: mc/3D/CitcomS/branches/compressible/lib/tracer_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/tracer_defs.h 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/tracer_defs.h 2007-01-19 20:14:41 UTC (rev 5838)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,10 +22,11 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
struct Tracer {
+
float *tracer_x;
float *tracer_y;
float *tracer_z;
@@ -33,5 +34,15 @@
float *x_space, *z_space, *y_space;
int NUM_TRACERS;
int LOCAL_NUM_TRACERS;
-int TEMP_TRACERS;
-}r;
+
+int *LOCAL_ELEMENT;
+
+float *THETA_LOC_ELEM;
+float *FI_LOC_ELEM;
+float *R_LOC_ELEM;
+
+float *THETA_LOC_ELEM_T;
+float *FI_LOC_ELEM_T;
+float *R_LOC_ELEM_T;
+
+} r;
Modified: mc/3D/CitcomS/branches/compressible/lib/viscosity_descriptions.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/viscosity_descriptions.h 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/lib/viscosity_descriptions.h 2007-01-19 20:14:41 UTC (rev 5838)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,18 +22,18 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
/* in this file define the contents of the VISC_OPT data structure
- which is used to store information used to create predefined
+ which is used to store information used to create predefined
viscosity fields, those determined from prior input, those
related to temperature/pressure/stress/anything else. */
struct VISC_OPT {
void (* update_viscosity)();
-
+
int update_allowed; /* determines whether visc field can evolve */
int EQUIVDD; /* Whatever the structure, average in the end */
int equivddopt;
@@ -41,13 +41,13 @@
int proflocy;
int SMOOTH;
int smooth_cycles;
-
+
char STRUCTURE[20]; /* which option to determine viscosity field, one of .... */
int FROM_SYSTEM;
int FROM_FILE;
int FROM_SPECS;
-
+
/* System ... */
int RHEOL; /* 1,2 */
int rheol_layers;
@@ -66,6 +66,14 @@
float freeze_thresh;
float freeze_value;
+ int channel;
+ int wedge;
+
+ float lv_min_radius;
+ float lv_max_radius;
+ float lv_channel_thickness;
+ float lv_reduction;
+
int MAX;
float max_value;
int MIN;
@@ -89,7 +97,7 @@
float weak_blobz[40];
float weak_blobwidth[40];
float weak_blobmag[40];
-
+
int weak_zones;
float weak_zonex1[40];
float weak_zoney1[40];
@@ -97,14 +105,14 @@
float weak_zonex2[40];
float weak_zoney2[40];
float weak_zonez2[40];
-
+
float weak_zonewidth[40];
float weak_zonemag[40];
-
+
int guess;
char old_file[100];
/* Specification info */
-
+
/* Prespecified viscosity parameters */
char VISC_OPT[20];
@@ -121,10 +129,10 @@
float cosx_epsilon;
float cosx_k;
int cosx_exp;
-
+
int EXPX;
float expx_epsilon;
-
+
/* MODULE BASED VISCOSITY VARIATIONS */
int RESDEPV;
@@ -133,5 +141,5 @@
int CHEMDEPV;
float CH0[40];
float CHEMeta0[40];
-
+
} viscosity;
Modified: mc/3D/CitcomS/branches/compressible/module/misc.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/module/misc.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/module/misc.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -281,7 +281,7 @@
E = (struct All_variables*)(PyCObject_AsVoidPtr(obj));
if(E->control.tracer==1)
- tracer_advection(E);
+ E->problem_tracer_advection(E);
Py_INCREF(Py_None);
return Py_None;
Modified: mc/3D/CitcomS/branches/compressible/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-01-19 19:44:49 UTC (rev 5837)
+++ mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-01-19 20:14:41 UTC (rev 5838)
@@ -265,6 +265,7 @@
num_perturb, fp);
}
else if (E->convection.tic_method == 1) {
+ getFloatProperty(properties, "half_space_age", E->convection.half_space_age, fp);
}
else if (E->convection.tic_method == 2) {
getFloatProperty(properties, "half_space_age", E->convection.half_space_age, fp);
@@ -694,6 +695,14 @@
getFloatVectorProperty(properties, "sdepv_expt",
E->viscosity.sdepv_expt, num_mat, fp);
+ getIntProperty(properties, "low_visc_channel", E->viscosity.channel, fp);
+ getIntProperty(properties, "low_visc_wedge", E->viscosity.wedge, fp);
+
+ getFloatProperty(properties, "lv_min_radius", E->viscosity.lv_min_radius, fp);
+ getFloatProperty(properties, "lv_max_radius", E->viscosity.lv_max_radius, fp);
+ getFloatProperty(properties, "lv_channel_thickness", E->viscosity.lv_channel_thickness, fp);
+ getFloatProperty(properties, "lv_reduction", E->viscosity.lv_reduction, fp);
+
getIntProperty(properties, "VMIN", E->viscosity.MIN, fp);
getFloatProperty(properties, "visc_min", E->viscosity.min_value, fp);
More information about the cig-commits
mailing list