[cig-commits] r14149 - in seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D: . citcoms_isotropic_no_crust
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Feb 25 10:31:51 PST 2009
Author: tan2
Date: 2009-02-25 10:31:51 -0800 (Wed, 25 Feb 2009)
New Revision: 14149
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/1D
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/build.mk
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_isotropic_no_crust.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_parm.h
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/crust
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/description.txt
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c
Removed:
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_prem_isotropic_no_crust/
Log:
Read (Vp, Vs, rho) from citcoms output as an earth model.
Note that no crust and no 1D reference model is used.
The interpolation part is unfinised yet.
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/1D
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/1D (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/1D 2009-02-25 18:31:51 UTC (rev 14149)
@@ -0,0 +1 @@
+link ../../1D_ref/blank/
\ No newline at end of file
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/1D
___________________________________________________________________
Name: svn:special
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/build.mk
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/build.mk (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/build.mk 2009-02-25 18:31:51 UTC (rev 14149)
@@ -0,0 +1,12 @@
+
+mantle_model_OBJECTS = \
+ $O/citcoms_isotropic_no_crust.o \
+ $O/read_citcoms_data.o \
+ $(EMPTY_MACRO)
+
+$O/citcoms_isotropic_no_crust.o: constants.h citcoms_isotropic_no_crust.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/citcoms_isotropic_no_crust.o ${FCFLAGS_f90} $S/MODELS/3D_mantle/citcoms_isotropic_no_crust/citcoms_isotropic_no_crust.f90
+
+$O/read_citcoms_data.o: citcoms_parm.h read_citcoms_data.c
+ ${CC} -c -o $O/read_citcoms_data.o ${CPPFLAGS} ${CFLAGS} $S/MODELS/3D_mantle/citcoms_isotropic_no_crust/read_citcoms_data.c
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_isotropic_no_crust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_isotropic_no_crust.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_isotropic_no_crust.f90 2009-02-25 18:31:51 UTC (rev 14149)
@@ -0,0 +1,294 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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 Softwait Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+
+!---------------------------
+
+subroutine get_model_properties(HONOR_1D_SPHERICAL_MOHO,ONE_CRUST, &
+ TRANSVERSE_ISOTROPY, &
+ ISOTROPIC_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ CRUSTAL,CASE_3D, &
+ ATTENUATION_3D)
+
+ implicit none
+
+ logical HONOR_1D_SPHERICAL_MOHO,ONE_CRUST,&
+ TRANSVERSE_ISOTROPY,&
+ ISOTROPIC_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,&
+ CRUSTAL,CASE_3D,&
+ ATTENUATION_3D
+
+ HONOR_1D_SPHERICAL_MOHO = .false.
+ ONE_CRUST = .true.
+
+ TRANSVERSE_ISOTROPY = .true.
+
+ ISOTROPIC_3D_MANTLE = .true.
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+
+ CRUSTAL = .false.
+ CASE_3D = .false.
+
+ ATTENUATION_3D = .false.
+
+end subroutine get_model_properties
+
+!---------------------------
+
+subroutine read_3d_mantle_model()
+
+ implicit none
+
+ ! coordinates of the four corner nodes in the current slice
+ double precision, dimension(4) :: x, y, z
+
+ call find_domain_corners(x, y, z)
+
+ call read_citcoms_mantle_model(x, y, z)
+
+end subroutine read_3d_mantle_model
+
+!---------------------------
+
+subroutine iso_mantle_model(radius, theta, phi, &
+ vpv, vph, vsv, vsh, rho, eta_aniso)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision, intent(in) :: radius, theta, phi
+ double precision, intent(out) :: vpv, vph, vsv, vsh, rho, eta_aniso
+ double precision inv_time
+
+ ! a static variable holding the result of previous search
+ integer, save :: elem
+ data elem /-1/
+
+ ! get (Vp, Vs, rho) in km/s and g/cm^3 from citcoms data
+ call get_velo_from_citcoms(radius, theta, phi, vpv, vsv, rho, elem)
+
+
+ ! non-dimensionalize
+ ! time scaling (s^{-1}) is done with inv_time
+ inv_time = dsqrt(PI * GRAV * RHOAV)
+
+ rho = rho * 1000.0d0 / RHOAV
+ vpv = vpv * 1000.0d0 / (R_EARTH * inv_time)
+ vph = vph * 1000.0d0 / (R_EARTH * inv_time)
+
+ vph = vpv
+ vsh = vsv
+
+ eta_aniso = one
+
+
+end subroutine iso_mantle_model
+
+!---------------------------
+
+subroutine aniso_mantle_model(r,theta,phi,rho, &
+ c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+
+ implicit none
+
+ double precision r,theta,phi
+ double precision rho
+ double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+ ! no-op
+
+end subroutine aniso_mantle_model
+
+!----------------------------------
+
+subroutine add_moho_topography(myrank,xelm,yelm,zelm,RMOHO,R220)
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision RMOHO,R220
+
+ ! no-np
+
+end subroutine add_moho_topography
+
+!----------------------------------
+
+subroutine add_topography_410_650(myrank,xelm,yelm,zelm,R220,R400,R670,R771)
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision R220,R400,R670,R771
+
+ ! no-op
+
+end subroutine add_topography_410_650
+
+!----------------------------------
+
+subroutine find_domain_corners(x, y, z)
+
+ implicit none
+
+ ! coordinates of the four corner nodes in the current slice
+ double precision, dimension(4) :: x, y, z
+
+ integer ichunk, iproc_xi, iproc_eta
+
+ !---------------------------------------------------------
+
+ call find_domain_partition(ichunk, iproc_xi, iproc_eta)
+
+ call find_corner_coord(ichunk, iproc_xi, iproc_eta, x, y, z)
+
+end subroutine find_domain_corners
+
+!----------------------------------
+
+subroutine find_domain_partition(ichunk, iproc_xi, iproc_eta)
+
+ implicit none
+
+ include "mpif.h"
+
+ integer, intent(out) :: ichunk, iproc_xi, iproc_eta
+
+ integer i, j, k
+ integer iprocnum, myrank, ier
+ integer NPROC_XI, NPROC_ETA, NCHUNKS
+ integer NPROC
+
+ ! the variables in the common block are for read only
+ common /for_citcoms1/ NPROC_XI, NPROC_ETA, NCHUNKS
+
+ !---------------------------------------------------------
+
+ ! find out the domain docomposition scheme
+ call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ier)
+
+ NPROC = NPROC_XI * NPROC_ETA
+
+ ! this do-loop is adapted from meshfem3D.f90
+ do i = 1,NCHUNKS
+ do j = 0,NPROC_ETA-1
+ do k = 0,NPROC_XI-1
+ iprocnum = (i-1)*NPROC + j * NPROC_XI + k
+ if(myrank == iprocnum) then
+ ichunk = i
+ iproc_eta = j
+ iproc_xi = k
+ ! write(*,*) iprocnum,ichunk,iproc_xi,iproc_eta
+ endif
+ enddo
+ enddo
+ enddo
+
+end subroutine find_domain_partition
+
+!----------------------------------
+
+subroutine find_corner_coord(ichunk, iproc_xi, iproc_eta, x, y, z)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: ichunk, iproc_xi, iproc_eta
+ double precision, dimension(4) :: x, y, z
+
+ ! local variables (some come from common blocks)
+ double precision, dimension(NGNOD) :: xelm, yelm, zelm, offset_x, offset_y, offset_z
+ double precision ANGULAR_WIDTH_XI_RAD, ANGULAR_WIDTH_ETA_RAD
+ double precision, dimension(NDIM,NDIM) :: rotation_matrix
+ double precision r_top, r_bottom, r_dummy
+
+ integer i, ilayer, nmesh_layers
+ integer NPROC_XI, NPROC_ETA, NCHUNKS
+
+ ! topology of the elements
+ integer, dimension(NGNOD) :: iaddx, iaddy, iaddz
+
+ ! the variables in the common block are for read only
+ common /for_citcoms1/ NPROC_XI, NPROC_ETA, NCHUNKS
+ common /for_citcoms2/ ANGULAR_WIDTH_XI_RAD, ANGULAR_WIDTH_ETA_RAD
+ common /for_citcoms3/ rotation_matrix
+
+ !
+ ! create a big element near the surface that covers the entire slice surface
+ !
+
+ ! create the shape of the corner nodes of a regular mesh element
+ call hex_nodes(iaddx, iaddy, iaddz)
+ offset_x(:) = iaddx(:) * HALF
+ offset_y(:) = iaddy(:) * HALF
+ offset_z(:) = iaddz(:) * HALF
+
+ ! r_top should be the raidus of the earth, r_bottom is unimportant
+ r_top = R_EARTH
+ r_bottom = 0.9d0 * R_EARTH
+
+ ! the actual layer is unimportant, as long as ilayer is not 1 (crust)
+ ilayer = 2
+ nmesh_layers = 10
+
+ !write(*,*) iproc_xi, iproc_eta, NPROC_XI, NPROC_ETA, ichunk
+
+ ! compute the transformed coordinates
+ call compute_coord_main_mesh(offset_x, offset_y, offset_z, &
+ xelm, yelm, zelm, &
+ ANGULAR_WIDTH_XI_RAD, ANGULAR_WIDTH_ETA_RAD, &
+ iproc_xi, iproc_eta, &
+ NPROC_XI, NPROC_ETA, 1, 1, &
+ r_top, r_bottom, &
+ 1, ilayer, ichunk, &
+ rotation_matrix, &
+ NCHUNKS, .false., nmesh_layers)
+
+ ! we only need the coordinate of the top nodes
+ x(:) = xelm(5:8)
+ y(:) = yelm(5:8)
+ z(:) = zelm(5:8)
+
+end subroutine find_corner_coord
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_parm.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_parm.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/citcoms_parm.h 2009-02-25 18:31:51 UTC (rev 14149)
@@ -0,0 +1,36 @@
+/*=====================================================================
+
+ S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+ --------------------------------------------------
+
+ Main authors: Dimitri Komatitsch and Jeroen Tromp
+ Seismological Laboratory, California Institute of Technology, USA
+ and University of Pau / CNRS / INRIA, France
+ (c) California Institute of Technology and University of Pau / CNRS / INRIA
+ February 2008
+
+ 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.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+=====================================================================*/
+
+
+/* filename (including full path) of citcoms domain output */
+const char citcoms_bound_filename[] = "/home/tan2/dv/trunk/fulltest.domain";
+
+/* filename base of citcoms output */
+const char citcoms_model_filename_base[] = "/home/tan2/dv/trunk/fulltest";
+
+/* time step of citcoms output */
+const int citcoms_step = 0;
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/crust
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/crust (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/crust 2009-02-25 18:31:51 UTC (rev 14149)
@@ -0,0 +1 @@
+link ../../crustal/none/
\ No newline at end of file
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/crust
___________________________________________________________________
Name: svn:special
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/description.txt
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/description.txt (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/description.txt 2009-02-25 18:31:51 UTC (rev 14149)
@@ -0,0 +1 @@
+The code reads the (theta, phi, r, rho, Vp, Vs) output from CitcomS, then interpolates to the GLL points.
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/3D/citcoms_isotropic_no_crust/read_citcoms_data.c 2009-02-25 18:31:51 UTC (rev 14149)
@@ -0,0 +1,1060 @@
+
+#include "config.h"
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <limits.h>
+#include <float.h>
+#include "mpi.h"
+#include "citcoms_parm.h"
+
+
+struct citcoms_cap {
+ double x[4];
+ double y[4];
+ double z[4];
+ double cost[4]; /* cos(theta) of this node */
+ double sint[4]; /* sin(theta) of this node */
+ double cosf[4]; /* cos(phi) of this node */
+ double sinf[4]; /* sin(phi) of this node */
+
+ double midx, midy, midz;
+ double disk_radius;
+};
+typedef struct citcoms_cap cap_t;
+
+
+struct citcoms_node {
+ double x;
+ double y;
+ double z;
+
+ double rho;
+ double vp;
+ double vs;
+};
+typedef struct citcoms_node node_t;
+
+
+struct citcoms_model {
+ double *rsqr; /* 1d radius^2 profile */
+
+ double *cost; /* cos(theta) of surface nodes */
+ double *sint; /* sin(theta) of surface nodes */
+ double *cosf; /* cos(phi) of surface nodes */
+ double *sinf; /* sin(phi) of surface nodes */
+
+ node_t *nodes;
+};
+typedef struct citcoms_model model_t;
+
+
+struct citcoms_domain {
+ int rank;
+
+ double rmin;
+ double rmax;
+
+ cap_t cap;
+
+ model_t model;
+
+ struct citcoms_domain *next;
+};
+typedef struct citcoms_domain domain_t;
+
+
+/* model static variable declaration */
+domain_t *first;
+
+/* citcoms local mesh geometry */
+int nox = 0;
+int noy = 0;
+int noz = 0;
+
+
+
+/************************************************************************
+ * Private functions first. Public functions at the end of the file.
+ ************************************************************************/
+
+
+/*
+ * Scan input str to get a long int vector *values. The vector length is from
+ * input len. The input str contains white-space seperated numbers.
+ */
+static int scan_long_vector(const char *str, int len, long *values)
+{
+ char *nptr, *endptr;
+ int i;
+
+ /* cast to avoid compiler warning */
+ nptr = endptr = (char *) str;
+
+ for (i = 0; i < len; ++i) {
+ values[i] = strtol(nptr, &endptr, 10);
+ if (nptr == endptr) {
+ /* error: no conversion is performed */
+ return -1;
+ }
+ nptr = endptr;
+ }
+
+ /** debug **
+ for (i = 0; i < len; ++i) fprintf(stderr, "%ld, ", values[i]);
+ fprintf(stderr, "\n");
+ /**/
+ return len;
+}
+
+
+/*
+ * From input file, read a line, which contains white-space seperated numbers
+ * of length num_columns, store the numbers in a long int array.
+ * Return num_columns on success, 0 on file reading error, and -1 on conversion error.
+ */
+static int read_long_vector(FILE *in, int num_columns, long *fields)
+{
+ char buffer[256], *p;
+
+ p = fgets(buffer, 255, in);
+ if (!p) {
+ return 0;
+ }
+
+ return scan_long_vector(buffer, num_columns, fields);
+}
+
+
+/*
+ * Scan input str to get a double vector *values. The vector length is from
+ * input len. The input str contains white-space seperated numbers.
+ */
+static int scan_double_vector(const char *str, int len, double *values)
+{
+ char *nptr, *endptr;
+ int i;
+
+ /* cast to avoid compiler warning */
+ nptr = endptr = (char *) str;
+
+ for (i = 0; i < len; ++i) {
+ values[i] = strtod(nptr, &endptr);
+ if (nptr == endptr) {
+ /* error: no conversion is performed */
+ return -1;
+ }
+ nptr = endptr;
+ }
+
+ /** debug **
+ for (i = 0; i < len; ++i) fprintf(stderr, "%e, ", values[i]);
+ fprintf(stderr, "\n");
+ /**/
+ return len;
+}
+
+
+/*
+ * From input file, read a line, which contains white-space seperated numbers
+ * of length num_columns, store the numbers in a double array.
+ * Return num_columns on success, 0 on file reading error, and -1 on conversion error.
+ */
+static int read_double_vector(FILE *in, int num_columns, double *fields)
+{
+ char buffer[256], *p;
+
+ p = fgets(buffer, 255, in);
+ if (!p) {
+ return 0;
+ }
+
+ return scan_double_vector(buffer, num_columns, fields);
+}
+
+
+static void create_model(model_t *model, const double *coord, const double *data)
+{
+ const double *x, *y, *z;
+ const double *rho, *vp, *vs;
+ int i, j, n;
+ const int nsf = nox * noy;
+ const int nno = nox * noy * noz;
+
+ /*** allocate space ****/
+ model->rsqr = malloc(noz *sizeof(double));
+
+ model->cost = malloc(nsf*sizeof(double));
+ model->sint = malloc(nsf*sizeof(double));
+ model->cosf = malloc(nsf*sizeof(double));
+ model->sinf = malloc(nsf*sizeof(double));
+
+ model->nodes = malloc(nno*sizeof(node_t));
+
+ if(model->rsqr == NULL ||
+ model->cost == NULL ||
+ model->sint == NULL ||
+ model->cosf == NULL ||
+ model->sinf == NULL ||
+ model->nodes == NULL) {
+ fprintf(stderr, "error while allocating memory '%s' 111\n",
+ __FILE__);
+ abort();
+ }
+
+ /*** store data ***/
+
+ /* some aliases for convenience */
+ x = &coord[0];
+ y = &coord[nno];
+ z = &coord[2*nno];
+ rho = &data[0];
+ vp = &data[nno];
+ vs = &data[2*nno];
+
+ /* the ordering of citcoms nodes is z first, then x, then y */
+ for(j=0; j<noz; j++) {
+ model->rsqr[j] = x[j]*x[j] + y[j]*y[j] + z[j]*z[j];
+ }
+
+ /* put node-related data in a struct to gain cache performance */
+ for(i=0, n=0; i<nsf; i++, n+=noz) {
+ double a, b;
+
+ a = sqrt(x[n]*x[n] + y[n]*y[n]);
+ b = sqrt(x[n]*x[n] + y[n]*y[n] + z[n]*z[n]);
+
+ model->cost[i] = z[n] / b;
+ model->sint[i] = a / b;
+ model->cosf[i] = x[n] / a;
+ model->sinf[i] = y[n] / a;
+
+ for(j=n; j<n+noz; j++) {
+ model->nodes[j].x = x[j];
+ model->nodes[j].y = y[j];
+ model->nodes[j].z = z[j];
+
+ model->nodes[j].rho = rho[j];
+ model->nodes[j].vp = vp[j];
+ model->nodes[j].vs = vs[j];
+ }
+ }
+
+ return;
+}
+
+
+static void new_cap_corners(cap_t *cap, const double x[4], const double y[4],
+ const double z[4])
+{
+ double r;
+ int i;
+
+ /* coordinates of corners */
+ memcpy(cap->x, x, 4*sizeof(double));
+ memcpy(cap->y, y, 4*sizeof(double));
+ memcpy(cap->z, z, 4*sizeof(double));
+
+ for(i=0; i<4; i++) {
+ double rad, rho;
+ rho = sqrt(x[i] * x[i] + y[i] * y[i]);
+ rad = sqrt(x[i] * x[i] + y[i] * y[i] * z[i] * z[i]);
+ cap->cost[i] = z[i] / rad;
+ cap->sint[i] = rho / rad;
+ cap->cosf[i] = x[i] / rho;
+ cap->sinf[i] = y[i] / rho;
+ }
+
+ /* coordinates of mid points, project to the unit sphere */
+ cap->midx = cap->midy = cap->midz = 0;
+ for(i=0; i<4; i++) {
+ cap->midx += x[i] * 0.25;
+ cap->midy += y[i] * 0.25;
+ cap->midz += z[i] * 0.25;
+ }
+ r = cap->midx * cap->midx + cap->midy * cap->midy + cap->midz * cap->midz;
+ cap->midx /= r;
+ cap->midy /= r;
+ cap->midz /= r;
+
+ /* radius of covering disk */
+ cap->disk_radius = 0;
+ for(i=0; i<4; i++) {
+ double tmp;
+
+ tmp = sqrt( (cap->midx - x[i]) * (cap->midx - x[i]) +
+ (cap->midy - y[i]) * (cap->midy - y[i]) +
+ (cap->midz - z[i]) * (cap->midz - z[i]) );
+ cap->disk_radius = (cap->disk_radius > tmp) ? cap->disk_radius : tmp;
+ }
+ return;
+}
+
+
+/* make a vector pointing from (x0,y0,z0) to (xf,yf,zf) */
+
+static void makevector(double vec[3], double xf, double yf, double zf,
+ double x0, double y0, double z0)
+{
+ vec[0] = xf - x0;
+ vec[1] = yf - y0;
+ vec[2] = zf - z0;
+ return;
+}
+
+
+/* cross product of vector A and B */
+
+static void crossit(double cross[3], const double A[3], const double B[3])
+{
+ cross[0] = A[1] * B[2] - A[2] * B[1];
+ cross[1] = A[2] * B[0] - A[0] * B[2];
+ cross[2] = A[0] * B[1] - A[1] * B[0];
+ return;
+}
+
+
+/* finds the radial component of a Cartesian vector */
+
+static double findradial(const double vec[3],
+ double cost, double sint,
+ double cosf, double sinf)
+{
+ double radialparti, radialpartj, radialpartk;
+ double radial;
+
+ radialparti = vec[0] * sint * cosf;
+ radialpartj = vec[1] * sint * sinf;
+ radialpartk = vec[2] * cost;
+
+ radial = radialparti + radialpartj + radialpartk;
+
+ return radial;
+}
+
+
+/***** ICHECK BOUNDS ******************************/
+/* */
+/* This function check if a test_point is bounded */
+/* by 4 corner nodes of a cap */
+/* This is done by: */
+/* 1) generate vectors from node to node */
+/* 2) generate vectors from each node to point */
+/* in question */
+/* 3) for each node, take cross product of vector */
+/* pointing to it from previous node and */
+/* vector from node to point in question */
+/* 4) Find radial components of all the cross */
+/* products. */
+/* 5) If all radial components are positive, */
+/* point is bounded by the 4 nodes */
+/* 6) If some radial components are close to 0 */
+/* point is on a boundary - adjust it an */
+/* epsilon amount for this analysis only */
+/* which will force it to lie in one element */
+/* or cap */
+/* */
+/* The 4 corner nodes must be ordered */
+/* counterclockwisely. */
+
+static int icheck_bounds(const double corner_x[4], const double corner_y[4],
+ const double corner_z[4],
+ const double corner_cost[4], const double corner_sint[4],
+ const double corner_cosf[4], const double corner_sinf[4],
+ double x, double y, double z)
+{
+ int ntries;
+ int i, icheck;
+
+ double v[4][3];
+ double rad[4];
+
+ double theta, phi;
+
+ const double eps = 1e-6;
+ const double tiny = 1e-15;
+
+ /* make vectors from corner node to corner node */
+ makevector(v[0], corner_x[0], corner_y[0], corner_z[0],
+ corner_x[3], corner_y[3], corner_z[3]);
+ for(i=0; i<3; i++)
+ makevector(v[i+1], corner_x[i+1], corner_y[i+1], corner_z[i+1],
+ corner_x[i], corner_y[i], corner_z[i]);
+
+ for(ntries=0; ntries<3; ntries++) {
+
+ for(i=0; i<4; i++) {
+ double p[3];
+ double cross[3];
+
+ /* make vectors from test point to node */
+ makevector(p, x, y, z, corner_x[i], corner_y[i], corner_z[i]);
+
+ /* Calculate cross products */
+ crossit(cross, v[i], p);
+
+ /* Calculate radial component of cross products */
+ rad[i] = findradial(cross, corner_cost[i], corner_sint[i],
+ corner_cosf[i], corner_sinf[i]);
+ }
+
+ /*
+ fprintf(stderr,"Rads: %f %f %f %f\n", rad[0], rad[1], rad[2], rad[3]);
+ fprintf(stderr,"Test Point: %f %f %f %f %f\n", x, y, z,
+ atan2(sqrt(x*x+y*y),z)*180/M_PI, atan2(y,x)*180/M_PI);
+ fprintf(stderr,"Corner points: 0: %f %f %f\n",
+ corner_x[0], corner_y[0], corner_z[0]);
+ fprintf(stderr,"Corner points: 1: %f %f %f\n",
+ corner_x[1], corner_y[1], corner_z[1]);
+ fprintf(stderr,"Corner points: 2: %f %f %f\n",
+ corner_x[2], corner_y[2], corner_z[2]);
+ fprintf(stderr,"Corner points: 3: %f %f %f\n",
+ corner_x[3], corner_y[3], corner_z[3]);
+ */
+
+ /* If all radial components are positive, point is bounded by the 4 nodes */
+ if (rad[0] >= 0 &&
+ rad[1] >= 0 &&
+ rad[2] >= 0 &&
+ rad[3] >= 0)
+ return 1;
+
+ /* If one radial component is significantly negative, point is not bounded */
+ if (rad[0] < -tiny ||
+ rad[1] < -tiny ||
+ rad[2] < -tiny ||
+ rad[3] < -tiny)
+ return 0;
+
+ /* Check if any radial components is zero (along a boundary), adjust if so */
+ /* Hopefully, this doesn't happen often, may be expensive */
+
+ theta = atan2(sqrt(x*x + y*y), z);
+ phi = atan2(y, x);
+
+ /* move toward the equator */
+ if (theta<=M_PI/2.0) {
+ theta += eps;
+ }
+ else {
+ theta -= eps;
+ }
+
+ phi += eps;
+
+ x = sin(theta) * cos(phi);
+ y = sin(theta) * sin(phi);
+ z = cos(theta);
+ }
+
+ /* even after nudging the point, it is still on the boundary, fail loudly */
+ fprintf(stderr,"Error(icheck_bounds) - too many tries\n");
+ fprintf(stderr,"Rads: %f %f %f %f\n", rad[0], rad[1], rad[2], rad[3]);
+ fprintf(stderr,"Test Point: %f %f %f %f %f\n", x, y, z,
+ theta*180/M_PI, phi*180/M_PI);
+ fprintf(stderr,"Corner points: 0: %f %f %f\n",
+ corner_x[0], corner_y[0], corner_z[0]);
+ fprintf(stderr,"Corner points: 1: %f %f %f\n",
+ corner_x[1], corner_y[1], corner_z[1]);
+ fprintf(stderr,"Corner points: 2: %f %f %f\n",
+ corner_x[2], corner_y[2], corner_z[2]);
+ fprintf(stderr,"Corner points: 3: %f %f %f\n",
+ corner_x[3], corner_y[3], corner_z[3]);
+ exit(10);
+
+ return 0;
+}
+
+
+/*
+ * In: one cap, consists of four points on a unit sphere and a convering disk; and
+ * one point on a unit sphere.
+ * Out: whether the point is inside the cap or not
+ */
+static int is_inside(const cap_t *cap, const double x, const double y, const double z)
+{
+ double dist;
+
+ dist = sqrt( (cap->midx - x) * (cap->midx - x) +
+ (cap->midy - y) * (cap->midy - y) +
+ (cap->midz - z) * (cap->midz - z) );
+
+ if(dist > cap->disk_radius)
+ return 0;
+
+ if(icheck_bounds(cap->x, cap->y, cap->z,
+ cap->cost, cap->sint, cap->cosf, cap->sinf,
+ x, y, z))
+ return 1;
+
+ return 0;
+}
+
+
+/*
+ * In: two caps, each consists of four points on a unit sphere and a convering disk
+ * Out: whether these caps overlap or not
+ */
+static int is_overlapping(const cap_t *cap0, const cap_t *cap1)
+{
+ double dist;
+ int i;
+
+ /* distance between cap centers */
+ dist = sqrt( (cap0->midx - cap1->midx) * (cap0->midx - cap1->midx) +
+ (cap0->midy - cap1->midy) * (cap0->midy - cap1->midy) +
+ (cap0->midz - cap1->midz) * (cap0->midz - cap1->midz) );
+
+ /* the centers are too far away then overlapping is impossible */
+ if(dist > (cap0->disk_radius + cap1->disk_radius))
+ return 0;
+
+ /* the center of one disk falls on another disk */
+ if(dist < cap0->disk_radius && dist < cap1->disk_radius)
+ return 1;
+
+ /* check each corner, is it inside the other cap? */
+ for(i=0; i<4; i++)
+ if(is_inside(cap0, cap1->x[i], cap1->y[i], cap1->z[i]))
+ return 1;
+
+ for(i=0; i<4; i++)
+ if(is_inside(cap1, cap0->x[i], cap0->y[i], cap0->z[i]))
+ return 1;
+
+ /* check each center, is it inside the other cap? */
+ if(icheck_bounds(cap0->x, cap0->y, cap0->z,
+ cap0->cost, cap0->sint, cap0->cosf, cap0->sinf,
+ cap1->midx, cap1->midy, cap1->midz))
+ return 1;
+
+ if(icheck_bounds(cap1->x, cap1->y, cap1->z,
+ cap1->cost, cap1->sint, cap1->cosf, cap1->sinf,
+ cap0->midx, cap0->midy, cap0->midz))
+ return 1;
+
+ return 0;
+}
+
+/* given an element number, return the node numbers of the bottom 4 nodes */
+static void bottom_nodes_of_element(int elem, int corners[4])
+{
+ int nx, ny, nz;
+ int lower_left_corner;
+ int i;
+
+ const int elxz = (nox - 1) * (noz - 1);
+
+ /* node connectivity ordered counterclockwisely */
+ const int conn[4] = {0, noz, noz*(nox+1), noz*nox};
+
+ /* The ordering scheme in citcoms is:
+ * elem = nz + nx*(noz-1) + ny*(nox-1)*(noz-1)
+ * node = nz + nx*noz + ny*nox*noz
+ */
+
+ nz = elem % (noz - 1);
+ ny = elem / elxz;
+ nx = (elem - ny * elxz) / (noz - 1);
+
+ /* the first node of this element */
+ lower_left_corner = elem + nx + ny * (nox + noz - 1);
+
+ for(i=0; i<4; i++)
+ corners[i] = lower_left_corner + conn[i];
+
+ return;
+}
+
+
+/* whether point (x, y, z) is inside element elem */
+static int is_inside_element(const model_t *model,
+ double x, double y, double z, int elem)
+{
+ double corner_x[4], corner_y[4], corner_z[4];
+ double corner_cost[4], corner_sint[4];
+ double corner_cosf[4], corner_sinf[4];
+ int i, j, corners[4];
+
+ bottom_nodes_of_element(elem, corners);
+
+ for(i=0; i<4; i++) {
+ corner_x[i] = model->nodes[corners[i]].x;
+ corner_y[i] = model->nodes[corners[i]].y;
+ corner_z[i] = model->nodes[corners[i]].z;
+
+ j = corners[i] / noz;
+ corner_cost[i] = model->cost[j];
+ corner_sint[i] = model->sint[j];
+ corner_cosf[i] = model->cosf[j];
+ corner_sinf[i] = model->sinf[j];
+ }
+
+ return icheck_bounds(corner_x, corner_y, corner_z,
+ corner_cost, corner_sint,
+ corner_cosf, corner_sinf,
+ x, y, z);
+}
+
+
+static int find_layer(const model_t *model,
+ double x, double y, double z)
+{
+ int i;
+ double r2 = x*x + y*y + z*z;
+
+ /* if point is below or above model domain */
+ if(r2 < model->rsqr[0] ||
+ r2 > model->rsqr[noz-1])
+ return -1;
+
+ /* return value: 0 to (noz-1) */
+ for(i=1; i<noz; i++) {
+ if(r2 < model->rsqr[i])
+ return i - 1;
+ }
+
+ return -1;
+}
+
+
+static int find_nearest_node(const model_t *model,
+ double x, double y, double z,
+ int nz)
+{
+ double dist2, min_dist2;
+ int i, n;
+ const int nno = nox * noy * noz;
+
+ min_dist2 = DBL_MAX;
+
+ /* loop over all nodes at this layer*/
+ n = -1;
+ for(i=nz; i<nno; i+=noz) {
+ node_t *node = &model->nodes[i];
+ dist2 = ( (node->x - x) * (node->x - x) +
+ (node->y - y) * (node->y - y) +
+ (node->z - z) * (node->z - z) );
+ if(dist2 < min_dist2) {
+ min_dist2 = dist2;
+ n = i;
+ }
+ }
+
+ return n;
+}
+
+
+static int ijk2elem(int nx, int ny, int nz)
+{
+ return nz + nx * (noz - 1) + ny * (nox - 1) * (noz - 1);
+}
+
+
+static int search_neighboring_elements(const model_t *model,
+ double x, double y, double z,
+ int nz, int nearest_node)
+{
+ int neighbor_elements[4];
+ int nx, ny;
+ int i;
+
+ const int nel = (nox-1) * (noy-1) * (noz-1);
+
+ ny = nearest_node / (nox * noz);
+ nx = (nearest_node - ny * nox * noz) / noz;
+
+ neighbor_elements[0] = ijk2elem(nx, ny, nz);
+ neighbor_elements[1] = ijk2elem(nx, ny-1, nz);
+ neighbor_elements[2] = ijk2elem(nx-1, ny-1, nz);
+ neighbor_elements[3] = ijk2elem(nx-1, ny, nz);
+
+ for(i=0; i<4; i++) {
+ if(neighbor_elements[i] < 0 || neighbor_elements[i] >= nel)
+ /* out of range */
+ continue;
+
+ if(is_inside_element(model, x, y, z, neighbor_elements[i]))
+ /* found */
+ return neighbor_elements[i];
+ }
+
+ /* not found */
+ return -1;
+}
+
+
+static int search_all_elements_in_layer(const model_t *model,
+ double x, double y, double z,
+ int nz)
+{
+ int nx, ny;
+ int elem;
+
+ for(nx=0; nx<nox-1; nx++)
+ for(ny=0; ny<noy-1; ny++) {
+ elem = ijk2elem(nx, ny, nz);
+ if(is_inside_element(model, x, y, z, elem))
+ /* found */
+ return elem;
+ }
+
+ /* not found */
+ return -1;
+}
+
+
+static int find_element(const model_t *model,
+ double x, double y, double z, int hint_elem)
+{
+ int nz, nearest_node, elem;
+
+ if(hint_elem > 0) {
+ /* check the hint_elem first */
+ if(is_inside_element(model, x, y, z, hint_elem))
+ /* found */
+ return hint_elem;
+ }
+
+ /* find corresponding radial layer */
+ nz = find_layer(model, x, y, z);
+ if(nz < 0)
+ /* not found */
+ return -1;
+
+ /* in this layer, find the nearest node to (x,y,z) */
+ nearest_node = find_nearest_node(model, x, y, z, nz);
+
+ /* search neighboring elements of the nearest node */
+ elem = search_neighboring_elements(model, x, y, z, nz, nearest_node);
+
+ if(elem >= 0)
+ /* found */
+ return elem;
+
+ /* a more expensive search */
+ elem = search_all_elements_in_layer(model, x, y, z, nz);
+
+ if(elem >= 0)
+ /* found */
+ return elem;
+
+
+ /* not found */
+ return -1;
+}
+
+
+/* convert (radius, theta, phi) on a unit sphere to (x, y, z) */
+static void rtp2xyz(double radius, double theta, double phi,
+ double *x, double *y, double *z)
+{
+ *x = radius * sin(theta) * cos(phi);
+ *y = radius * sin(theta) * sin(phi);
+ *z = radius * cos(theta);
+ return;
+}
+
+
+static void get_shape_functions(const model_t *model, int elem,
+ double x, double y, double z,
+ double shp[8])
+{
+ //XXX
+ return;
+}
+
+
+static void interpolate_velocities(const model_t *model, const double shp[8],
+ double *vp, double *vs, double *rho)
+{
+ int i;
+ for(i=0; i<8; i++) {
+ //XXX: do something
+ }
+ return;
+}
+
+
+static void get_velocities(const model_t *model, int elem,
+ double x, double y, double z,
+ double *vp, double *vs, double *rho)
+{
+ double shp[8];
+
+ get_shape_functions(model, elem, x, y, z, shp);
+ interpolate_velocities(model, shp, vp, vs, rho);
+
+ return;
+}
+
+
+
+
+/************************************************************************
+ * Public functions
+ ************************************************************************/
+
+void FC_FUNC_(read_citcoms_mantle_model, READ_CITCOMS_MANTLE_MODEL)
+ (const double x[4], const double y[4], const double z[4])
+{
+ const int ncolumns = 11;
+ int myrank;
+ int i, buffer_size;
+ int total_citcoms_proc;
+ double *buffer;
+ double cx[4], cy[4], cz[4];
+ FILE *fp;
+ cap_t my_cap;
+ domain_t *curr, *prev;
+
+ MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+ new_cap_corners(&my_cap, x, y, z);
+
+
+ /***********************************************************************
+ * Read citcoms domain bound. Rank-0 processor reads the file and
+ * broadcast to other processors.
+ * See citcoms manual for file format
+ ***********************************************************************/
+
+ if(myrank == 0) {
+ long ltmp[2];
+
+ fp = fopen(citcoms_bound_filename, "r");
+ if(fp == NULL) {
+ fprintf(stderr, "error while opening file '%s'\n",
+ citcoms_bound_filename);
+ abort();
+ }
+
+ /* read header line */
+ if(read_long_vector(fp, 2, ltmp) != 2) {
+ fprintf(stderr, "error while reading first line of '%s'\n",
+ citcoms_bound_filename);
+ abort();
+ }
+
+ if(ltmp[1] != ncolumns) {
+ fprintf(stderr, "'%s' contains wrong number of columns\n",
+ citcoms_bound_filename);
+ abort();
+ }
+
+ total_citcoms_proc = (int) ltmp[0];
+ }
+
+ MPI_Bcast(&total_citcoms_proc, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ buffer_size = total_citcoms_proc * ncolumns;
+ buffer = malloc(buffer_size * sizeof(double));
+ if(buffer == NULL) {
+ fprintf(stderr, "error while allocating memory '%s' 1\n",
+ __FILE__);
+ abort();
+ }
+
+ if(myrank == 0) {
+ /* read content of domain bound */
+ for(i=0; i<total_citcoms_proc; i++) {
+ if(read_double_vector(fp, ncolumns, &(buffer[i*ncolumns])) != ncolumns) {
+ fprintf(stderr, "error while reading %d-th line of '%s'\n",
+ i, citcoms_bound_filename);
+ abort();
+ }
+ }
+ fclose(fp);
+ }
+
+ MPI_Bcast(buffer, buffer_size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+
+
+
+ /***********************************************************************
+ * Find out which citcoms domain overlaps with my domain.
+ ***********************************************************************/
+
+ first = malloc(sizeof(domain_t));
+ curr = first;
+ prev = NULL;
+
+ for(i=0; i<total_citcoms_proc; i++) {
+ double *tmp = &(buffer[i*ncolumns]);
+ double ctheta, cphi;
+ cap_t citcoms_cap;
+ int j;
+
+ for(j=0; j<4; j++) {
+ ctheta = tmp[1 + 2*j];
+ cphi = tmp[2 + 2*j];
+ rtp2xyz(1.0, ctheta, cphi, &(cx[j]), &(cy[j]), &(cz[j]));
+ }
+
+ new_cap_corners(&citcoms_cap, cx, cy, cz);
+
+ if(is_overlapping(&my_cap, &citcoms_cap)) {
+ curr->rank = (int) tmp[0];
+ curr->rmin = tmp[9];
+ curr->rmax = tmp[10];
+ memcpy(&(curr->cap), &citcoms_cap, sizeof(cap_t));
+ curr->next = malloc(sizeof(domain_t));
+ prev = curr;
+ curr = curr->next;
+ }
+ }
+
+ if(prev == NULL) {
+ fprintf(stderr, "rank=%d: cannot find overllaping domain '%s' 2\n",
+ myrank, __FILE__);
+ abort();
+ }
+
+ /* this pointer is over-allocated */
+ free(curr);
+ prev->next = NULL;
+
+ free(buffer);
+
+ /* count how many overlapping domains */
+ curr = first;
+ i = 0;
+ do {
+ i++;
+ } while (curr = curr->next);
+ fprintf(stderr, "rank=%d, overlaps with %d citcoms domains\n", myrank, i);
+ MPI_Barrier(MPI_COMM_WORLD);
+
+
+ /***********************************************************************
+ * Read the mantle model from those citcoms domains
+ ***********************************************************************/
+
+ curr = first;
+ do {
+ const char fmt1[] = "%s.coord_bin.%d";
+ const char fmt2[] = "%s.seismic.%d.%d";
+ char filename1[256];
+ char filename2[256];
+ FILE *f1, *f2;
+ int header[4];
+ int nno, size;
+ double *coord, *seismic;
+
+ snprintf(filename1, 255, fmt1, curr->rank);
+ f1 = fopen(filename1, "rb");
+ if(f1 == NULL) {
+ fprintf(stderr, "error while opening file '%s'\n",
+ filename1);
+ abort();
+ }
+
+ snprintf(filename2, 255, fmt2, curr->rank, citcoms_step);
+ f2 = fopen(filename2, "rb");
+ if(f2 == NULL) {
+ fprintf(stderr, "error while opening file '%s'\n",
+ filename2);
+ abort();
+ }
+
+ /* read header */
+ fread(header, sizeof(int), 4, f1);
+ if(header[3] != INT_MAX) {
+ fprintf(stderr, "error while reading file header '%s'\n",
+ filename1);
+ abort();
+ }
+ nox = header[0];
+ noy = header[1];
+ noz = header[2];
+ nno = nox * noy * noz;
+
+ coord = malloc(nno*3*sizeof(double));
+ seismic = malloc(nno*3*sizeof(double));
+ if(coord == NULL || seismic == NULL) {
+ fprintf(stderr, "error while allocating memory rank=%d\n", myrank);
+ abort();
+ }
+
+ /* read coordinates */
+ size = fread(coord, sizeof(double), nno*3, f1);
+ if(size != nno*3) {
+ fprintf(stderr, "error while reading file '%s'\n",
+ filename1);
+ abort();
+ }
+
+ /* read (rho, vp, vs) */
+ size = fread(seismic, sizeof(double), nno*3, f2);
+ if(size != nno*3) {
+ fprintf(stderr, "error while reading file '%s'\n",
+ filename1);
+ abort();
+ }
+
+ /* store data */
+ create_model(&(curr->model), coord, seismic);
+
+ free(seismic);
+ free(coord);
+
+ fclose(f2);
+ fclose(f1);
+ } while (curr = curr->next);
+
+ fprintf(stderr, "rank=%d, 10\n", myrank);
+ MPI_Barrier(MPI_COMM_WORLD);
+
+ return;
+}
+
+
+void FC_FUNC_(get_velo_from_citcoms, GET_VELO_FROM_CITCOMS)
+ (const double *radius, const double *theta, const double *phi,
+ double *vpv, double *vsv, double *rho,
+ int *hint_elem)
+{
+ int elem;
+ double x, y, z;
+ domain_t *curr;
+
+ curr = first;
+ do {
+
+ if(*radius < curr->rmin || *radius > curr->rmax)
+ /* not in this citcoms domain */
+ continue;
+
+ rtp2xyz(*radius, *theta, *phi, &x, &y, &z);
+
+ if(! is_inside(&(curr->cap), x, y, z))
+ /* not in this citcoms domain */
+ continue;
+
+ elem = find_element(&(curr->model), x, y, z, *hint_elem);
+
+ if(elem < 0)
+ /* search failed!!! */
+ break;
+
+ /* search is successful */
+ get_velocities(&(curr->model), elem, x, y, z,
+ vpv, vsv, rho);
+
+ /* the next search is very likely to be in the same element,
+ * save the current result as the hint of next search */
+ *hint_elem = elem;
+
+ return;
+
+ } while (curr = curr->next);
+
+ /* failed search */
+ fprintf(stderr, "Cannot find corresponding citcoms element! (%.15e, %.15e, %.15e)\n",
+ *radius, *theta, *phi);
+ abort();
+
+ return;
+}
+
+
More information about the CIG-COMMITS
mailing list