https://github.com/geodynamics/citcoms
Tip revision: 304266d8447efe5ee990624f615b250c633079ef authored by Leif Strand on 14 February 2006, 23:27:17 UTC
Merged r2605, "Updated NEWS and ChangeLog".
Merged r2605, "Updated NEWS and ChangeLog".
Tip revision: 304266d
Tracer_advection.c
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*
*<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 for CitcomS
*/
#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,nn;
int m;
float idummy,xdummy,ydummy,zdummy;
FILE *fp, *fp1;
int nox,noy,noz,gnox,gnoy,gnoz;
char aaa[100];
char output_file[255];
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);
}
fp1=fopen("coorr.dat","r");
if (fp1 == NULL) {
fprintf(E->fp,"(Tracer_advection #2) Cannot open %s\n","coorr.dat");
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));
fscanf(fp1,"%s",aaa);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
for(i=1;i<=gnox;i++)
{
/*
E->Tracer.x_space[i]=E->sphere.cap[m].theta[1]+(E->sphere.cap[m].theta[3]-E->sphere.cap[m].theta[1])/(gnox-1)*(i-1);
*/
fscanf(fp1,"%d %f",&nn,&(E->Tracer.x_space[i]));
if(E->parallel.me==0) printf("%s %f\n","E->Tracer.x_space[i]", E->Tracer.x_space[i]);
}
fscanf(fp1,"%s",aaa);
for(k=1;k<=gnoy;k++)
{
/*
E->Tracer.y_space[k]=E->sphere.cap[m].fi[1]+(E->sphere.cap[m].fi[3]-E->sphere.cap[m].fi[1])/(gnoy-1)*(k-1);
*/
fscanf(fp1,"%d %f",&nn,&(E->Tracer.y_space[k]));
/*
if(E->parallel.me==0) printf("%s %f\n","E->Tracer.y_space[k]", E->Tracer.y_space[k]);
*/
}
fscanf(fp1,"%s",aaa);
for(j=1;j<=gnoz;j++)
{
/*
E->Tracer.z_space[j]=E->sphere.ri+(E->sphere.ro-E->sphere.ri)/(gnoz-1)*(j-1);
*/
fscanf(fp1,"%d %f",&nn,&(E->Tracer.z_space[j]));
/*
if(E->parallel.me==0) printf("%s %f\n","E->Tracer.z_space[j]", E->Tracer.z_space[j]);
*/
}
} /* end of m */
for(i=1;i<=E->Tracer.NUM_TRACERS;i++)
{
fscanf(fp,"%f %f %f %f",
&idummy, &xdummy, &ydummy, &zdummy);
E->Tracer.itcolor[i]=idummy; E->Tracer.tracer_x[i]=xdummy;
E->Tracer.tracer_y[i]=ydummy;E->Tracer.tracer_z[i]=zdummy;
}
return;
}
void tracer_advection(E)
struct All_variables *E;
{
float find_age_in_MY();
int n,loc;
int i,j,k,m,ii,jj,kk;
int n_x,n_y,n_z;
int node1,node2,node3,node4,node5,node6,node7,node8;
int nno,nox,noy,noz,gnox,gnoy,gnoz,nodeg,nodel;
int nxs,nys,nzs,proc,nprocxl,nprocyl,nproczl;
int iteration;
int root;
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[1][4];
float time;
MPI_Comm world;
float xmin,xmax,ymin,ymax,zmin,zmax;
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;
/*** comment by Tan2 1/25/2005
This part of code find the bouding box of the global mesh.
Not checked.
***/
for(m=1;m<=E->sphere.caps_per_proc;m++) {
if(E->sphere.cap[m].theta[1]<=E->sphere.cap[m].theta[3]) {
xmin=E->sphere.cap[m].theta[1];
xmax=E->sphere.cap[m].theta[3];
}
else {
xmin=E->sphere.cap[m].theta[3];
xmax=E->sphere.cap[m].theta[1];
}
if(E->sphere.cap[m].fi[1]<=E->sphere.cap[m].fi[3]) {
ymin=E->sphere.cap[m].fi[1];
ymax=E->sphere.cap[m].fi[3];
}
else {
ymin=E->sphere.cap[m].fi[3];
ymax=E->sphere.cap[m].fi[1];
}
zmin=E->sphere.ri;
zmax=E->sphere.ro;
} /* end of m */
/*** comment by Tan2 1/25/2005
Copy the velocity array. Is it necessary?
***/
for(m=1;m<=E->sphere.caps_per_proc;m++) {
for(i=0;i<=nno-1;i++)
for(j=1;j<=3;j++) {
E->V[m][j][i]=E->sphere.cap[m].V[j][i+1];
}
}
//time=find_age_in_MY(E); /* NOT TESTED CPC 6/25/00 */
root=0;
/*** comment by Tan2 1/25/2005
Gather the global velocity to root processor from all processors.
Seems the code is designed for multi-processor run, although the
tracer advection only runs in root processor.
***/
for(m=1;m<=E->sphere.caps_per_proc;m++) {
MPI_Gather(E->V[m][1], nno, MPI_FLOAT, E->GV1[m][1], nno, MPI_FLOAT, root, world);
MPI_Gather(E->V[m][2], nno, MPI_FLOAT, E->GV1[m][2], nno, MPI_FLOAT, root, world);
MPI_Gather(E->V[m][3], nno, MPI_FLOAT, E->GV1[m][3], nno, MPI_FLOAT, root, world);
}
loc=0;
nprocxl=E->parallel.nprocx;
nprocyl=E->parallel.nprocy;
nproczl=E->parallel.nprocz;
/*** comment by Tan2 1/25/2005
Combine the gathered velocity to get the global velocity field.
***/
if(E->parallel.me==root) {
for(m=1;m<=E->sphere.caps_per_proc;m++)
for (k=1;k<=nprocyl;k++)
for (j=1;j<=nprocxl;j++)
for (i=1;i<=nproczl;i++) {
proc = i-1 + (j-1)*nproczl+(k-1)*nprocxl*nproczl;
nxs = (j-1)*(nox-1) + 1;
nys = (k-1)*(noy-1) + 1;
nzs = (i-1)*(noz-1) + 1;
for (kk=1;kk<=noy;kk++)
for (jj=1;jj<=nox;jj++)
for (ii=1;ii<=noz;ii++) {
nodeg = ii + nzs -1 + gnoz*
(nxs+jj-2)
+ (nys+kk-2)*gnox*gnoz;
nodel = (kk-1)*nox*noz + (jj-1)*noz + ii;
E->GV[m][1][nodeg] = E->GV1[m][1][nodel+loc-1];
E->GV[m][2][nodeg] = E->GV1[m][2][nodel+loc-1];
E->GV[m][3][nodeg] = E->GV1[m][3][nodel+loc-1];
}
loc += nno;
}
}
/*** comment by Tan2 1/25/2005
Now, root processor has velocity field of the whole mesh and starts
to advect the tracers.
***/
if(E->parallel.me==root) {
for(n=1;n<=E->Tracer.NUM_TRACERS;n++) {
n_x=0;
n_y=0;
n_z=0;
/*** comment by Tan2 1/25/2005
itcolor is like the starting age of the tracer. Disable it temperory.
***/
//if(E->Tracer.itcolor[n]>=time) {
if(1) {
/* 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<gnox;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<gnoz;j++) {
if(z_tmp >= E->Tracer.z_space[j] && z_tmp <= E->Tracer.z_space[j+1]) {
tr_dz=z_tmp-E->Tracer.z_space[j];
dz=E->Tracer.z_space[j+1]-E->Tracer.z_space[j];
n_z=j;
}
}
for(k=1;k<gnoy;k++) {
if(y_tmp >= E->Tracer.y_space[k] && y_tmp <= E->Tracer.y_space[k+1]) {
tr_dy=y_tmp-E->Tracer.y_space[k];
dy=E->Tracer.y_space[k+1]-E->Tracer.y_space[k];
n_y=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)*gnoz + (n_y-1)*gnoz*gnox;
node2 = n_z + n_x*gnoz + (n_y-1)*gnoz*gnox;
node3 = n_z+1 + (n_x-1)*gnoz + (n_y-1)*gnoz*gnox;
node4 = n_z+1 + n_x*gnoz + (n_y-1)*gnoz*gnox;
node5 = n_z + (n_x-1)*gnoz + n_y*gnoz*gnox;
node6 = n_z + n_x*gnoz + n_y*gnoz*gnox;
node7 = n_z+1 + (n_x-1)*gnoz + n_y*gnoz*gnox;
node8 = n_z+1 + n_x*gnoz + n_y*gnoz*gnox;
/*** 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 */
} /* if E->Tracer.itcolor */
} /* end of tracer loop */
}
/* make sure that the new position of the tracers within the box */
for(n=1;n<=E->Tracer.NUM_TRACERS;n++) {
if(E->Tracer.tracer_x[n]>xmax) E->Tracer.tracer_x[n]=2.0*xmax-E->Tracer.tracer_x[n];
if(E->Tracer.tracer_x[n]<xmin) E->Tracer.tracer_x[n]=2.0*xmin-E->Tracer.tracer_x[n];
if(E->Tracer.tracer_y[n]>ymax) E->Tracer.tracer_y[n]=2.0*ymax-E->Tracer.tracer_y[n];
if(E->Tracer.tracer_y[n]<ymin) E->Tracer.tracer_y[n]=2.0*ymin-E->Tracer.tracer_y[n];
if(E->Tracer.tracer_z[n]>zmax) E->Tracer.tracer_z[n]=2.0*zmax-E->Tracer.tracer_z[n];
if(E->Tracer.tracer_z[n]<zmin) E->Tracer.tracer_z[n]=2.0*zmin-E->Tracer.tracer_z[n];
}
return;
}