https://github.com/geodynamics/citcoms
Revision df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC, committed by Thorsten Becker on 08 December 2010, 18:43:58 UTC
GMT/ggrd compile before).
1 parent cf7495e
Tip revision: df821c06952128812763ec5b7a944cd178ca5ed3 authored by Thorsten Becker on 08 December 2010, 18:43:58 UTC
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Added definition of TRUE (1) and FALSE (0) in case undefined (was defined for
Tip revision: df821c0
Regional_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>
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#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"
#include "composition_related.h"
#include "parallel_related.h"
static void write_trace_instructions(struct All_variables *E);
static void make_mesh_ijk(struct All_variables *E);
static void put_lost_tracers(struct All_variables *E,
int *send_size, double *send,
int kk, int j);
static void put_found_tracers(struct All_variables *E,
int recv_size, double *recv,
int j);
int isearch_neighbors(double *array, int nsize,
double a, int hint);
int isearch_all(double *array, int nsize, double a);
void regional_tracer_setup(struct All_variables *E)
{
char output_file[255];
void get_neighboring_caps();
double CPU_time0();
double begin_time = CPU_time0();
/* Some error control */
if (E->sphere.caps_per_proc>1) {
fprintf(stderr,"This code does not work for multiple caps per processor!\n");
parallel_process_termination();
}
/* open tracing output file */
sprintf(output_file,"%s.tracer_log.%d",E->control.data_file,E->parallel.me);
E->trace.fpt=fopen(output_file,"w");
/* reset statistical counters */
E->trace.istat_isend=0;
E->trace.istat_iempty=0;
E->trace.istat_elements_checked=0;
E->trace.istat1=0;
/* some obscure initial parameters */
/* This parameter specifies how close a tracer can get to the boundary */
E->trace.box_cushion=0.00001;
/* Determine number of tracer quantities */
/* advection_quantites - those needed for advection */
E->trace.number_of_basic_quantities=12;
/* extra_quantities - used for flavors, composition, etc. */
/* (can be increased for additional science i.e. tracing chemistry */
E->trace.number_of_extra_quantities = 0;
if (E->trace.nflavors > 0)
E->trace.number_of_extra_quantities += 1;
E->trace.number_of_tracer_quantities =
E->trace.number_of_basic_quantities +
E->trace.number_of_extra_quantities;
/* Fixed positions in tracer array */
/* Flavor is always in extraq position 0 */
/* Current coordinates are always kept in basicq positions 0-5 */
/* Other positions may be used depending on science being done */
/* Some error control regarding size of pointer arrays */
if (E->trace.number_of_basic_quantities>99) {
fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of basic in tracer_defs.h\n");
fflush(E->trace.fpt);
parallel_process_termination();
}
if (E->trace.number_of_extra_quantities>99) {
fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of extraq in tracer_defs.h\n");
fflush(E->trace.fpt);
parallel_process_termination();
}
if (E->trace.number_of_tracer_quantities>99) {
fprintf(E->trace.fpt,"ERROR(initialize_trace)-increase 2nd position size of rlater in tracer_defs.h\n");
fflush(E->trace.fpt);
parallel_process_termination();
}
write_trace_instructions(E);
/* The bounding box of neiboring processors */
get_neighboring_caps(E);
make_mesh_ijk(E);
if (E->composition.on)
composition_setup(E);
fprintf(E->trace.fpt, "Tracer intiailization takes %f seconds.\n",
CPU_time0() - begin_time);
return;
}
/**** WRITE TRACE INSTRUCTIONS ***************/
static void write_trace_instructions(struct All_variables *E)
{
int i;
fprintf(E->trace.fpt,"\nTracing Activated! (proc: %d)\n",E->parallel.me);
fprintf(E->trace.fpt," Allen K. McNamara 12-2003\n\n");
if (E->trace.ic_method==0) {
fprintf(E->trace.fpt,"Generating New Tracer Array\n");
fprintf(E->trace.fpt,"Tracers per element: %d\n",E->trace.itperel);
}
if (E->trace.ic_method==1) {
fprintf(E->trace.fpt,"Reading tracer file %s\n",E->trace.tracer_file);
}
if (E->trace.ic_method==2) {
fprintf(E->trace.fpt,"Read individual tracer files\n");
}
fprintf(E->trace.fpt,"Number of tracer flavors: %d\n", E->trace.nflavors);
if (E->trace.nflavors && E->trace.ic_method==0) {
fprintf(E->trace.fpt,"Initialized tracer flavors by: %d\n", E->trace.ic_method_for_flavors);
if (E->trace.ic_method_for_flavors == 0) {
fprintf(E->trace.fpt,"Layered tracer flavors\n");
for (i=0; i<E->trace.nflavors-1; i++)
fprintf(E->trace.fpt,"Interface Height: %d %f\n",i,E->trace.z_interface[i]);
}
#ifdef USE_GGRD
else if((E->trace.ic_method_for_flavors == 1)||(E->trace.ic_method_for_flavors == 99)) {
/* ggrd modes 1 and 99 (99 is override for restart) */
fprintf(stderr,"ggrd regional flavors not implemented\n");
fprintf(E->trace.fpt,"ggrd not implemented et for regional, flavor method= %d\n",
E->trace.ic_method_for_flavors);
fflush(E->trace.fpt);
parallel_process_termination();
}
#endif
else {
fprintf(E->trace.fpt,"Sorry-This IC methods for Flavors are Unavailable %d\n",E->trace.ic_method_for_flavors);
fflush(E->trace.fpt);
parallel_process_termination();
}
}
for (i=0; i<E->trace.nflavors-2; i++) {
if (E->trace.z_interface[i] < E->trace.z_interface[i+1]) {
fprintf(E->trace.fpt,"Sorry - The %d-th z_interface is smaller than the next one.\n", i);
fflush(E->trace.fpt);
parallel_process_termination();
}
}
/* more obscure stuff */
fprintf(E->trace.fpt,"Box Cushion: %f\n",E->trace.box_cushion);
fprintf(E->trace.fpt,"Number of Basic Quantities: %d\n",
E->trace.number_of_basic_quantities);
fprintf(E->trace.fpt,"Number of Extra Quantities: %d\n",
E->trace.number_of_extra_quantities);
fprintf(E->trace.fpt,"Total Number of Tracer Quantities: %d\n",
E->trace.number_of_tracer_quantities);
if (E->trace.itracer_warnings==0) {
fprintf(E->trace.fpt,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
fprintf(stderr,"\n WARNING EXITS ARE TURNED OFF! TURN THEM ON!\n");
fflush(E->trace.fpt);
}
write_composition_instructions(E);
return;
}
static void make_mesh_ijk(struct All_variables *E)
{
int m,i,j,k,node;
int nox,noy,noz;
nox=E->lmesh.nox;
noy=E->lmesh.noy;
noz=E->lmesh.noz;
E->trace.x_space=(double*) malloc(nox*sizeof(double));
E->trace.y_space=(double*) malloc(noy*sizeof(double));
E->trace.z_space=(double*) malloc(noz*sizeof(double));
/***comment by Vlad 1/26/2005
reading the local mesh coordinate
***/
for(m=1;m<=E->sphere.caps_per_proc;m++) {
for(i=0;i<nox;i++)
E->trace.x_space[i]=E->sx[m][1][i*noz+1];
for(j=0;j<noy;j++)
E->trace.y_space[j]=E->sx[m][2][j*nox*noz+1];
for(k=0;k<noz;k++)
E->trace.z_space[k]=E->sx[m][3][k+1];
} /* end of m */
/* debug *
for(i=0;i<nox;i++)
fprintf(E->trace.fpt, "i=%d x=%e\n", i, E->trace.x_space[i]);
for(j=0;j<noy;j++)
fprintf(E->trace.fpt, "j=%d y=%e\n", j, E->trace.y_space[j]);
for(k=0;k<noz;k++)
fprintf(E->trace.fpt, "k=%d z=%e\n", k, E->trace.z_space[k]);
/**
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 0));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 1));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 2));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 3));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.7, 4));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.56, 0));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.56, 1));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.56, 2));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.99, 2));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.99, 3));
fprintf(stderr, "%d\n", isearch_neighbors(E->trace.z_space, noz, 0.99, 4));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.5));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 1.1));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.55));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 1.0));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.551));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.99));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.7));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.75));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.775));
fprintf(stderr, "%d\n", isearch_all(E->trace.z_space, noz, 0.7750001));
parallel_process_termination();
/**/
return;
}
/********** IGET ELEMENT *****************************************/
/* */
/* This function returns the the real element for a given point. */
/* Returns -99 in not in this cap. */
/* iprevious_element, if known, is the last known element. If */
/* it is not known, input a negative number. */
int regional_iget_element(struct All_variables *E,
int m, int iprevious_element,
double dummy1, double dummy2, double dummy3,
double theta, double phi, double rad)
{
int e, i, j, k;
int ii, jj, kk;
int elx, ely, elz;
elx = E->lmesh.elx;
ely = E->lmesh.ely;
elz = E->lmesh.elz;
//TODO: take care of south west bound
/* Search neighboring elements if the previous element is known */
if (iprevious_element > 0) {
e = iprevious_element - 1;
k = e % elz;
i = (e / elz) % elx;
j = e / (elz*elx);
ii = isearch_neighbors(E->trace.x_space, elx+1, theta, i);
jj = isearch_neighbors(E->trace.y_space, ely+1, phi, j);
kk = isearch_neighbors(E->trace.z_space, elz+1, rad, k);
if (ii>=0 && jj>=0 && kk>=0)
return jj*elx*elz + ii*elz + kk + 1;
}
/* Search all elements if either the previous element is unknown */
/* or failed to find in the neighboring elements */
ii = isearch_all(E->trace.x_space, elx+1, theta);
jj = isearch_all(E->trace.y_space, ely+1, phi);
kk = isearch_all(E->trace.z_space, elz+1, rad);
if (ii<0 || jj<0 || kk<0)
return -99;
return jj*elx*elz + ii*elz + kk + 1;
}
/* array is an ordered array of length nsize */
/* return an index i, such that array[i] <= a < array[i+1] */
/* return -1 if not found. */
/* Note that -1 is returned if a == array[nsize-1] */
int isearch_all(double *array, int nsize, double a)
{
int high, i, low;
/* check the min/max bound */
if ((a < array[0]) || (a >= array[nsize-1]))
return -1;
/* binary search */
for (low=0, high=nsize-1; high-low>1;) {
i = (high+low) / 2;
if ( a < array[i] ) high = i;
else low = i;
}
return low;
}
/* Similar the isearch_all(), but with a hint */
int isearch_neighbors(double *array, int nsize,
double a, int hint)
{
/* search the nearest neighbors only */
const int number_of_neighbors = 3;
int neighbors[5];
int n, i;
neighbors[0] = hint;
neighbors[1] = hint-1;
neighbors[2] = hint+1;
neighbors[3] = hint-2;
neighbors[4] = hint+2;
/**/
for (n=0; n<number_of_neighbors; n++) {
i = neighbors[n];
if ((i >= 0) && (i < nsize-1) &&
(a >= array[i]) && (a < array[i+1]))
return i;
}
return -1;
}
/* */
/* This function serves to determine if a point lies within */
/* a given cap */
/* */
int regional_icheck_cap(struct All_variables *E, int icap,
double theta, double phi, double rad, double junk)
{
double theta_min, theta_max;
double phi_min, phi_max;
/* corner 2 is the north-west corner */
/* corner 4 is the south-east corner */
theta_min = E->trace.theta_cap[icap][2];
theta_max = E->trace.theta_cap[icap][4];
phi_min = E->trace.phi_cap[icap][2];
phi_max = E->trace.phi_cap[icap][4];
if ((theta >= theta_min) && (theta < theta_max) &&
(phi >= phi_min) && (phi < phi_max))
return 1;
//TODO: deal with south west bounds
return 0;
}
void regional_get_shape_functions(struct All_variables *E,
double shp[9], int nelem,
double theta, double phi, double rad)
{
int e, i, j, k;
int elx, ely, elz;
double tr_dx, tr_dy, tr_dz;
double dx, dy, dz;
double volume;
elx = E->lmesh.elx;
ely = E->lmesh.ely;
elz = E->lmesh.elz;
e = nelem - 1;
k = e % elz;
i = (e / elz) % elx;
j = e / (elz*elx);
/*** comment by Tan2 1/25/2005
Find the element that contains the tracer.
node(i) tracer node(i+1)
| * |
<----------->
tr_dx
<-------------------------------->
dx
***/
tr_dx = theta - E->trace.x_space[i];
dx = E->trace.x_space[i+1] - E->trace.x_space[i];
tr_dy = phi - E->trace.y_space[j];
dy = E->trace.y_space[j+1] - E->trace.y_space[j];
tr_dz = rad - E->trace.z_space[k];
dz = E->trace.z_space[k+1] - E->trace.z_space[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 */
volume = dx*dz*dy;
shp[1] = (dx-tr_dx) * (dy-tr_dy) * (dz-tr_dz) / volume;
shp[2] = tr_dx * (dy-tr_dy) * (dz-tr_dz) / volume;
shp[3] = tr_dx * tr_dy * (dz-tr_dz) / volume;
shp[4] = (dx-tr_dx) * tr_dy * (dz-tr_dz) / volume;
shp[5] = (dx-tr_dx) * (dy-tr_dy) * tr_dz / volume;
shp[6] = tr_dx * (dy-tr_dy) * tr_dz / volume;
shp[7] = tr_dx * tr_dy * tr_dz / volume;
shp[8] = (dx-tr_dx) * tr_dy * tr_dz / volume;
/** debug **
fprintf(E->trace.fpt, "dr=(%e,%e,%e) tr_dr=(%e,%e,%e)\n",
dx, dy, dz, tr_dx, tr_dy, tr_dz);
fprintf(E->trace.fpt, "shp: %e %e %e %e %e %e %e %e\n",
shp[1], shp[2], shp[3], shp[4], shp[5], shp[6], shp[7], shp[8]);
fprintf(E->trace.fpt, "sum(shp): %e\n",
shp[1]+ shp[2]+ shp[3]+ shp[4]+ shp[5]+ shp[6]+ shp[7]+ shp[8]);
fflush(E->trace.fpt);
/**/
return;
}
double regional_interpolate_data(struct All_variables *E,
double shp[9], double data[9])
{
int n;
double result = 0;
for(n=1; n<=8; n++)
result += data[n] * shp[n];
return result;
}
/******** GET VELOCITY ***************************************/
void regional_get_velocity(struct All_variables *E,
int m, int nelem,
double theta, double phi, double rad,
double *velocity_vector)
{
void velo_from_element_d();
double shp[9], VV[4][9], tmp;
int n, d, node;
const int sphere_key = 0;
/* get shape functions at (theta, phi, rad) */
regional_get_shape_functions(E, shp, nelem, theta, phi, rad);
/* get cartesian velocity */
velo_from_element_d(E, VV, m, nelem, sphere_key);
/*** comment by Tan2 1/25/2005
Interpolate the velocity on the tracer position
***/
for(d=1; d<=3; d++)
velocity_vector[d] = 0;
for(d=1; d<=3; d++) {
for(n=1; n<=8; n++)
velocity_vector[d] += VV[d][n] * shp[n];
}
/** debug **
for(d=1; d<=3; d++) {
fprintf(E->trace.fpt, "VV: %e %e %e %e %e %e %e %e: %e\n",
VV[d][1], VV[d][2], VV[d][3], VV[d][4],
VV[d][5], VV[d][6], VV[d][7], VV[d][8],
velocity_vector[d]);
}
tmp = 0;
for(n=1; n<=8; n++)
tmp += E->sx[m][1][E->ien[m][nelem].node[n]] * shp[n];
fprintf(E->trace.fpt, "THETA: %e -> %e\n", theta, tmp);
fflush(E->trace.fpt);
/**/
return;
}
void regional_keep_within_bounds(struct All_variables *E,
double *x, double *y, double *z,
double *theta, double *phi, double *rad)
{
void sphere_to_cart();
int changed = 0;
if (*theta > E->control.theta_max - E->trace.box_cushion) {
*theta = E->control.theta_max - E->trace.box_cushion;
changed = 1;
}
if (*theta < E->control.theta_min + E->trace.box_cushion) {
*theta = E->control.theta_min + E->trace.box_cushion;
changed = 1;
}
if (*phi > E->control.fi_max - E->trace.box_cushion) {
*phi = E->control.fi_max - E->trace.box_cushion;
changed = 1;
}
if (*phi < E->control.fi_min + E->trace.box_cushion) {
*phi = E->control.fi_min + E->trace.box_cushion;
changed = 1;
}
if (*rad > E->sphere.ro - E->trace.box_cushion) {
*rad = E->sphere.ro - E->trace.box_cushion;
changed = 1;
}
if (*rad < E->sphere.ri + E->trace.box_cushion) {
*rad = E->sphere.ri + E->trace.box_cushion;
changed = 1;
}
if (changed)
sphere_to_cart(E, *theta, *phi, *rad, x, y, z);
return;
}
void regional_lost_souls(struct All_variables *E)
{
/* This part only works if E->sphere.caps_per_proc==1 */
const int j = 1;
int lev = E->mesh.levmax;
int i, d, kk;
int max_send_size, isize, itemp_size;
int ngbr_rank[6+1];
double bounds[3][2];
double *send[2];
double *recv[2];
void expand_tracer_arrays();
int icheck_that_processor_shell();
int ipass;
MPI_Status status[4];
MPI_Request request[4];
double CPU_time0();
double begin_time = CPU_time0();
E->trace.istat_isend = E->trace.ilater[j];
/* the bounding box */
for (d=0; d<E->mesh.nsd; d++) {
bounds[d][0] = E->sx[j][d+1][1];
bounds[d][1] = E->sx[j][d+1][E->lmesh.nno];
}
/* set up ranks for neighboring procs */
/* if ngbr_rank is -1, there is no neighbor on this side */
ipass = 1;
for (kk=1; kk<=6; kk++) {
if (E->parallel.NUM_PASS[lev][j].bound[kk] == 1) {
ngbr_rank[kk] = E->parallel.PROCESSOR[lev][j].pass[ipass];
ipass++;
}
else {
ngbr_rank[kk] = -1;
}
}
/* debug *
for (kk=1; kk<=E->trace.istat_isend; kk++) {
fprintf(E->trace.fpt, "tracer#=%d xx=(%g,%g,%g)\n", kk,
E->trace.rlater[j][0][kk],
E->trace.rlater[j][1][kk],
E->trace.rlater[j][2][kk]);
}
for (d=0; d<E->mesh.nsd; d++) {
fprintf(E->trace.fpt, "bounds(dim=%d) = (%e, %e)\n",
d, bounds[d][0], bounds[d][1]);
}
for (kk=1; kk<=6; kk++) {
fprintf(E->trace.fpt, "pass=%d neighbor_rank=%d\n",
kk, ngbr_rank[kk]);
}
fflush(E->trace.fpt);
parallel_process_sync(E);
/**/
/* Allocate Maximum Memory to Send Arrays */
max_send_size = max(2*E->trace.ilater[j], E->trace.ntracers[j]/100);
itemp_size = max_send_size * E->trace.number_of_tracer_quantities;
if ((send[0] = (double *)malloc(itemp_size*sizeof(double)))
== NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (u388)\n");
fflush(E->trace.fpt);
exit(10);
}
if ((send[1] = (double *)malloc(itemp_size*sizeof(double)))
== NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (u389)\n");
fflush(E->trace.fpt);
exit(10);
}
for (d=0; d<E->mesh.nsd; d++) {
int original_size = E->trace.ilater[j];
int idb;
int kk = 1;
int isend[2], irecv[2];
isend[0] = isend[1] = 0;
/* move out-of-bound tracers to send array */
while (kk<=E->trace.ilater[j]) {
double coord;
/* Is the tracer within the bounds in the d-th dimension */
coord = E->trace.rlater[j][d][kk];
if (coord < bounds[d][0]) {
put_lost_tracers(E, &(isend[0]), send[0], kk, j);
}
else if (coord >= bounds[d][1]) {
put_lost_tracers(E, &(isend[1]), send[1], kk, j);
}
else {
/* check next tracer */
kk++;
}
/* reallocate send if size too small */
if ((isend[0] > max_send_size - 5) ||
(isend[1] > max_send_size - 5)) {
isize = max_send_size + max_send_size/4 + 10;
itemp_size = isize * E->trace.number_of_tracer_quantities;
if ((send[0] = (double *)realloc(send[0],
itemp_size*sizeof(double)))
== NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (s4)\n");
fflush(E->trace.fpt);
exit(10);
}
if ((send[1] = (double *)realloc(send[1],
itemp_size*sizeof(double)))
== NULL) {
fprintf(E->trace.fpt,"Error(lost souls)-no memory (s5)\n");
fflush(E->trace.fpt);
exit(10);
}
fprintf(E->trace.fpt,"Expanding physical memory of send to "
"%d from %d\n",
isize, max_send_size);
max_send_size = isize;
}
} /* end of while kk */
/* check the total # of tracers is conserved */
if ((isend[0] + isend[1] + E->trace.ilater[j]) != original_size) {
fprintf(E->trace.fpt, "original_size: %d, rlater_size: %d, "
"send_size: %d\n",
original_size, E->trace.ilater[j], kk);
}
/** debug **
for (i=0; i<2; i++) {
for (kk=0; kk<isend[i]; kk++) {
fprintf(E->trace.fpt, "dim:%d side:%d kk=%d coord[kk]=%e\n",
d, i, kk,
send[i][kk*E->trace.number_of_tracer_quantities+d]);
}
}
fflush(E->trace.fpt);
/**/
/* Send info to other processors regarding number of send tracers */
/* check whether there is a neighbor in this pass*/
idb = 0;
for (i=0; i<2; i++) {
int target_rank;
kk = d*2 + i + 1;
target_rank = ngbr_rank[kk];
if (target_rank >= 0) {
MPI_Isend(&isend[i], 1, MPI_INT, target_rank,
11, E->parallel.world, &request[idb++]);
MPI_Irecv(&irecv[i], 1, MPI_INT, target_rank,
11, E->parallel.world, &request[idb++]);
}
else {
irecv[i] = 0;
}
} /* end of for i */
/* Wait for non-blocking calls to complete */
MPI_Waitall(idb, request, status);
/** debug **
for (i=0; i<2; i++) {
int target_rank;
kk = d*2 + i + 1;
target_rank = ngbr_rank[kk];
if (target_rank >= 0) {
fprintf(E->trace.fpt, "%d: %d send %d to proc %d\n",
d, i, isend[i], target_rank);
fprintf(E->trace.fpt, "%d: %d recv %d from proc %d\n",
d, i, irecv[i], target_rank);
}
}
parallel_process_sync(E);
/**/
/* Allocate memory in receive arrays */
for (i=0; i<2; i++) {
isize = irecv[i] * E->trace.number_of_tracer_quantities;
itemp_size = max(1, isize);
if ((recv[i] = (double *)malloc(itemp_size*sizeof(double)))
== NULL) {
fprintf(E->trace.fpt, "Error(lost souls)-no memory (c721)\n");
fflush(E->trace.fpt);
exit(10);
}
}
/* Now, send the tracers to proper procs */
idb = 0;
for (i=0; i<2; i++) {
int target_rank;
kk = d*2 + i + 1;
target_rank = ngbr_rank[kk];
if (target_rank >= 0) {
isize = isend[i] * E->trace.number_of_tracer_quantities;
MPI_Isend(send[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
isize = irecv[i] * E->trace.number_of_tracer_quantities;
MPI_Irecv(recv[i], isize, MPI_DOUBLE, target_rank,
12, E->parallel.world, &request[idb++]);
}
}
/* Wait for non-blocking calls to complete */
MPI_Waitall(idb, request, status);
/** debug **
for (i=0; i<2; i++) {
for (kk=1; kk<=irecv[i]; kk++) {
fprintf(E->trace.fpt, "recv: %d %e %e %e\n",
kk,
recv[i][(kk-1)*E->trace.number_of_tracer_quantities],
recv[i][(kk-1)*E->trace.number_of_tracer_quantities+1],
recv[i][(kk-1)*E->trace.number_of_tracer_quantities+2]);
}
}
fflush(E->trace.fpt);
parallel_process_sync(E);
/**/
/* put the received tracers */
for (i=0; i<2; i++) {
put_found_tracers(E, irecv[i], recv[i], j);
}
free(recv[0]);
free(recv[1]);
} /* end of for d */
/* rlater should be empty by now */
if (E->trace.ilater[j] > 0) {
fprintf(E->trace.fpt, "Error(regional_lost_souls) lost tracers\n");
for (kk=1; kk<=E->trace.ilater[j]; kk++) {
fprintf(E->trace.fpt, "lost #%d xx=(%e, %e, %e)\n", kk,
E->trace.rlater[j][0][kk],
E->trace.rlater[j][1][kk],
E->trace.rlater[j][2][kk]);
}
fflush(E->trace.fpt);
exit(10);
}
/* Free Arrays */
free(send[0]);
free(send[1]);
E->trace.lost_souls_time += CPU_time0() - begin_time;
return;
}
static void put_lost_tracers(struct All_variables *E,
int *send_size, double *send,
int kk, int j)
{
int ilast_tracer, isend_position, ipos;
int pp;
/* move the tracer from rlater to send */
isend_position = (*send_size) * E->trace.number_of_tracer_quantities;
for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
ipos = isend_position + pp;
send[ipos] = E->trace.rlater[j][pp][kk];
}
(*send_size)++;
/* eject the tracer from rlater */
ilast_tracer = E->trace.ilater[j];
for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++) {
E->trace.rlater[j][pp][kk] = E->trace.rlater[j][pp][ilast_tracer];
}
E->trace.ilater[j]--;
return;
}
/****************************************************************/
/* Put the received tracers in basiq & extraq, if within bounds */
/* Otherwise, append to rlater for sending to another proc */
static void put_found_tracers(struct All_variables *E,
int recv_size, double *recv,
int j)
{
void expand_tracer_arrays();
void expand_later_array();
int icheck_processor_shell();
int kk, pp;
int ipos, ilast, inside, iel;
double theta, phi, rad;
for (kk=0; kk<recv_size; kk++) {
ipos = kk * E->trace.number_of_tracer_quantities;
theta = recv[ipos];
phi = recv[ipos + 1];
rad = recv[ipos + 2];
/* check whether this tracer is inside this proc */
/* check radius first, since it is cheaper */
inside = icheck_processor_shell(E, j, rad);
if (inside == 1)
inside = regional_icheck_cap(E, 0, theta, phi, rad, rad);
else
inside = 0;
/** debug **
fprintf(E->trace.fpt, "kk=%d, inside=%d, xx=(%e, %e, %e)\n",
kk, inside, theta, phi, rad);
fprintf(E->trace.fpt, "before: %d %d\n",
E->trace.ilater[j], E->trace.ntracers[j]);
/**/
if (inside) {
E->trace.ntracers[j]++;
ilast = E->trace.ntracers[j];
if (E->trace.ntracers[j] > (E->trace.max_ntracers[j]-5))
expand_tracer_arrays(E, j);
for (pp=0; pp<E->trace.number_of_basic_quantities; pp++)
E->trace.basicq[j][pp][ilast] = recv[ipos+pp];
ipos += E->trace.number_of_basic_quantities;
for (pp=0; pp<E->trace.number_of_extra_quantities; pp++)
E->trace.extraq[j][pp][ilast] = recv[ipos+pp];
/* found the element */
iel = regional_iget_element(E, j, -99, 0, 0, 0, theta, phi, rad);
if (iel<1) {
fprintf(E->trace.fpt, "Error(regional lost souls) - "
"element not here?\n");
fprintf(E->trace.fpt, "theta, phi, rad: %f %f %f\n",
theta, phi, rad);
fflush(E->trace.fpt);
exit(10);
}
E->trace.ielement[j][ilast] = iel;
}
else {
if (E->trace.ilatersize[j]==0) {
E->trace.ilatersize[j]=E->trace.max_ntracers[j]/5;
for (kk=0;kk<E->trace.number_of_tracer_quantities;kk++) {
if ((E->trace.rlater[j][kk]=(double *)malloc(E->trace.ilatersize[j]*sizeof(double)))==NULL) {
fprintf(E->trace.fpt,"AKM(put_found_tracers)-no memory (%d)\n",kk);
fflush(E->trace.fpt);
exit(10);
}
}
} /* end first particle initiating memory allocation */
E->trace.ilater[j]++;
ilast = E->trace.ilater[j];
if (E->trace.ilater[j] > (E->trace.ilatersize[j]-5))
expand_later_array(E, j);
for (pp=0; pp<E->trace.number_of_tracer_quantities; pp++)
E->trace.rlater[j][pp][ilast] = recv[ipos+pp];
} /* end of if-else */
/** debug **
fprintf(E->trace.fpt, "after: %d %d\n",
E->trace.ilater[j], E->trace.ntracers[j]);
fflush(E->trace.fpt);
/**/
} /* end of for kk */
return;
}
Computing file changes ...