[cig-commits] r9128 - cs/avm/trunk
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Jan 24 12:06:13 PST 2008
Author: tan2
Date: 2008-01-24 12:06:13 -0800 (Thu, 24 Jan 2008)
New Revision: 9128
Added:
cs/avm/trunk/avm-partition.c
cs/avm/trunk/avm.f90
cs/avm/trunk/write_CitcomS_tracers.f90
Log:
Initial import from Leif Strand.
* write_CitcomS_tracers.f90: reads binary SPECFEM mesher files (procXXXXXX_reg1_solver_data_2.bin); outputs a single tracer file with columns theta, phi, radius, rank, and iglob -- one line is generated for each GLL point in the crust and mantle (it ignores inner core and outer core).
* avm-partition.c: reads an arbitrary number of input files with the following format:
header (ignored)
rho kappa mu rank iglob
rho kappa mu rank iglob
rho kappa mu rank iglob
...
It partitions the data into NPROC files according to 'rank':
proc000000_reg1_avm.txt
proc000001_reg1_avm.txt
proc000002_reg1_avm.txt
...
Each such file has the following format:
iglob rho kappa mu
iglob rho kappa mu
iglob rho kappa mu
* avm.f90: reads each of the "procXXXXXX_reg1_avm.txt" files in sequence and saves rho, kappa, mu into the binary SPECFEM mesher files (procXXXXXX_reg1_solver_data_1.bin). I haven't tested this program; I couldn't think of an easy way to verify that it's working.
--Leif
Added: cs/avm/trunk/avm-partition.c
===================================================================
--- cs/avm/trunk/avm-partition.c (rev 0)
+++ cs/avm/trunk/avm-partition.c 2008-01-24 20:06:13 UTC (rev 9128)
@@ -0,0 +1,74 @@
+
+/* support files larger than 2GB on Linux */
+#define _FILE_OFFSET_BITS 64
+
+#include <stdio.h>
+#include <stdlib.h>
+
+
+#define IREGION_CRUST_MANTLE 1
+
+
+int main(int argc, char **argv)
+{
+ FILE **out; int nout;
+ int i;
+
+ if (argc < 2) {
+ fprintf(stderr, "usage: %s FILE...\n", argv[0]);
+ return 1;
+ }
+
+ nout = 0;
+ out = 0;
+ for (i = 1; i < argc; ++i) {
+ char *ifilename;
+ FILE *in;
+ int c;
+ double rho, kappa, mu, rank, glob;
+
+ ifilename = argv[i];
+ in = fopen(ifilename, "r");
+ if (!in) {
+ fprintf(stderr, "%s: fopen: `%s': ", argv[0], ifilename);
+ perror("");
+ return 1;
+ }
+
+ /* skip header (?) */
+ do {
+ c = fgetc(in);
+ if (c == EOF) {
+ fprintf(stderr, "%s: bad input file: `%s'\n", argv[0], ifilename);
+ return 1;
+ }
+ } while (c != '\n');
+
+ while (fscanf(in, "%lf %lf %lf %lf %lf\n", &rho, &kappa, &mu, &rank, &glob) == 5) {
+ int irank, nproc;
+
+ /* open more output files as needed */
+ irank = (int)rank;
+ nproc = irank + 1;
+ if (nproc > nout) {
+ out = (FILE **)realloc(out, nproc * sizeof(FILE *));
+ for ( ; nout < nproc; ++nout) {
+ char ofilename[32];
+ sprintf(ofilename, "proc%06d_reg%d_avm.txt", nout, IREGION_CRUST_MANTLE);
+ out[nout] = fopen(ofilename, "w");
+ }
+ }
+
+ fprintf(out[irank], "%d %lf %lf %lf\n", (int)glob, rho, kappa, mu);
+ }
+
+ fclose(in);
+ }
+
+ for (--nout; nout >= 0; --nout) {
+ fclose(out[nout]);
+ }
+ free(out);
+
+ return 0;
+}
Added: cs/avm/trunk/avm.f90
===================================================================
--- cs/avm/trunk/avm.f90 (rev 0)
+++ cs/avm/trunk/avm.f90 2008-01-24 20:06:13 UTC (rev 9128)
@@ -0,0 +1,244 @@
+!=====================================================================
+!
+! 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, France
+! (c) California Institute of Technology and University of Pau / CNRS, November 2007
+!
+! 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.
+!
+!=====================================================================
+
+program avm
+ ! read arrays created by the mesher
+
+ implicit none
+
+ include 'mpif.h'
+ include "constants.h"
+ include "precision.h"
+
+ ! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+ ! -------- arrays specific to each region here -----------
+
+ ! ----------------- crust, mantle and oceans ---------------------
+
+ ! mesh parameters
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: &
+ xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle
+
+ ! arrays for isotropic elements stored only where needed to save space
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
+ rhostore_crust_mantle,kappavstore_crust_mantle,muvstore_crust_mantle
+
+ ! arrays for anisotropic elements stored only where needed to save space
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_TISO_MANTLE) :: &
+ kappahstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle
+
+ ! arrays for full anisotropy only when needed
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: &
+ c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
+ c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle, &
+ c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle, &
+ c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle, &
+ c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
+ c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
+ c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle
+
+ ! local to global mapping
+ integer, dimension(NSPEC_CRUST_MANTLE) :: idoubling_crust_mantle
+
+ ! mass matrix
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rmass_crust_mantle
+
+ ! ------------------------------------
+
+ ! additional mass matrix for ocean load
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
+
+ ! Stacey
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STACEY) :: rho_vp_crust_mantle,rho_vs_crust_mantle
+
+ ! parameters read from parameter file
+ logical OCEANS, ABSORBING_CONDITIONS
+ character(len=150) LOCAL_PATH
+ ! processor identification
+ character(len=150) prname
+
+ ! flags to read kappa and mu and anisotropy arrays in regions where needed
+ logical READ_KAPPA_MU,READ_TISO
+ ! proc numbers for MPI
+ integer myrank
+
+ integer iglob
+ integer i,j,k
+
+ ! AVM files
+ integer ipoint,ispec
+ real(kind=CUSTOM_REAL) rho,kappa,mu
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: rho_avm,kappa_avm,mu_avm
+ logical mask_ibool(NGLOB_CRUST_MANTLE)
+
+ OCEANS = .false.
+ ABSORBING_CONDITIONS = .false.
+ LOCAL_PATH = 'scratch'
+
+
+ do myrank=0,NPROCTOT_VAL-1
+
+ ! create the name for the database of the current slide and region
+ call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
+
+ ! crust and mantle
+
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ READ_KAPPA_MU = .false.
+ READ_TISO = .false.
+ else
+ READ_KAPPA_MU = .true.
+ READ_TISO = .true.
+ endif
+
+ call read_arrays_solver(IREGION_CRUST_MANTLE,myrank, &
+ rho_vp_crust_mantle,rho_vs_crust_mantle, &
+ xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
+ rhostore_crust_mantle,kappavstore_crust_mantle,muvstore_crust_mantle, &
+ kappahstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle, &
+ NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE,NSPECMAX_ANISO_MANTLE, &
+ c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
+ c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle, &
+ c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle, &
+ c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle, &
+ c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
+ c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
+ c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
+ ibool_crust_mantle,idoubling_crust_mantle,rmass_crust_mantle,rmass_ocean_load, &
+ NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE, &
+ READ_KAPPA_MU,READ_TISO,TRANSVERSE_ISOTROPY_VAL,ANISOTROPIC_3D_MANTLE_VAL, &
+ ANISOTROPIC_INNER_CORE_VAL,OCEANS,LOCAL_PATH,ABSORBING_CONDITIONS)
+
+ ! -------------------------------------------
+ open(unit=IIN,file=trim(prname)//'avm.txt',status='old',action='read')
+
+ mask_ibool(:) = .true.
+
+ do ipoint=1,NGLOB_CRUST_MANTLE
+ read(IIN,*) iglob, rho, kappa, mu
+
+ if (iglob < 0 .or. iglob > NGLOB_CRUST_MANTLE) stop 'invalid global point number'
+ rho_avm(iglob) = rho
+ kappa_avm(iglob) = kappa
+ mu_avm(iglob) = mu
+ mask_ibool(iglob) = .false.
+ enddo
+
+ if (count(mask_ibool(:)) /= 0) stop 'uninitialized GLL points'
+
+ do ispec = 1, NSPEC_CRUST_MANTLE
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ rhostore_crust_mantle(i,j,k,ispec) = rho_avm(iglob)
+ kappavstore_crust_mantle(i,j,k,ispec) = kappa_avm(iglob)
+ muvstore_crust_mantle(i,j,k,ispec) = mu_avm(iglob)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ close(IIN)
+ ! -------------------------------------------
+
+ open(unit=27,file=trim(prname)//'solver_data_1.bin',status='unknown',form='unformatted')
+
+ write(27) xix_crust_mantle
+ write(27) xiy_crust_mantle
+ write(27) xiz_crust_mantle
+ write(27) etax_crust_mantle
+ write(27) etay_crust_mantle
+ write(27) etaz_crust_mantle
+ write(27) gammax_crust_mantle
+ write(27) gammay_crust_mantle
+ write(27) gammaz_crust_mantle
+
+ write(27) rhostore_crust_mantle
+ write(27) kappavstore_crust_mantle
+
+ if(READ_KAPPA_MU) write(27) muvstore_crust_mantle
+
+ if(TRANSVERSE_ISOTROPY_VAL .and. READ_TISO) then
+ write(27) kappahstore_crust_mantle
+ write(27) muhstore_crust_mantle
+ write(27) eta_anisostore_crust_mantle
+ endif
+
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ write(27) c11store_crust_mantle
+ write(27) c12store_crust_mantle
+ write(27) c13store_crust_mantle
+ write(27) c14store_crust_mantle
+ write(27) c15store_crust_mantle
+ write(27) c16store_crust_mantle
+ write(27) c22store_crust_mantle
+ write(27) c23store_crust_mantle
+ write(27) c24store_crust_mantle
+ write(27) c25store_crust_mantle
+ write(27) c26store_crust_mantle
+ write(27) c33store_crust_mantle
+ write(27) c34store_crust_mantle
+ write(27) c35store_crust_mantle
+ write(27) c36store_crust_mantle
+ write(27) c44store_crust_mantle
+ write(27) c45store_crust_mantle
+ write(27) c46store_crust_mantle
+ write(27) c55store_crust_mantle
+ write(27) c56store_crust_mantle
+ write(27) c66store_crust_mantle
+ endif
+
+ ! Stacey
+ if(ABSORBING_CONDITIONS) then
+
+ write(27) rho_vp_crust_mantle
+ write(27) rho_vs_crust_mantle
+
+ endif
+
+ ! mass matrix
+ write(27) rmass_crust_mantle
+
+ ! additional ocean load mass matrix if oceans and if we are in the crust
+ if(OCEANS) write(27) rmass_ocean_load
+
+ close(27)
+
+ enddo
+
+end program avm
Added: cs/avm/trunk/write_CitcomS_tracers.f90
===================================================================
--- cs/avm/trunk/write_CitcomS_tracers.f90 (rev 0)
+++ cs/avm/trunk/write_CitcomS_tracers.f90 2008-01-24 20:06:13 UTC (rev 9128)
@@ -0,0 +1,87 @@
+!=====================================================================
+!
+! 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, France
+! (c) California Institute of Technology and University of Pau / CNRS, November 2007
+!
+! 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.
+!
+!=====================================================================
+
+program write_CitcomS_tracers
+
+ implicit none
+
+ include 'mpif.h'
+ include "constants.h"
+ include "precision.h"
+
+ ! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: &
+ xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle
+ character(len=150) LOCAL_PATH, prname
+ integer myrank, iglob
+
+ ! tracers
+ integer num_tracers
+ real(kind=CUSTOM_REAL) theta,phi,radius
+ integer flavor1,flavor2
+
+
+ LOCAL_PATH = 'scratch'
+
+ open(unit=IOUT,file='tracer.dat',status='unknown')
+
+ num_tracers = NPROCTOT_VAL * NGLOB_CRUST_MANTLE
+ write(IOUT,*) num_tracers
+
+ do myrank=0,NPROCTOT_VAL-1
+
+ ! create the name for the database of the current slice
+ call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
+
+ ! read coordinates of the mesh
+ open(unit=IIN,file=trim(prname)//'solver_data_2.bin',status='old',action='read',form='unformatted')
+ read(IIN) xstore_crust_mantle
+ read(IIN) ystore_crust_mantle
+ read(IIN) zstore_crust_mantle
+ close(IIN)
+
+ do iglob=1,NGLOB_CRUST_MANTLE
+ ! convert to spherical coordinates
+ call xyz_2_rthetaphi( &
+ xstore_crust_mantle(iglob), &
+ ystore_crust_mantle(iglob), &
+ zstore_crust_mantle(iglob), &
+ radius, &
+ theta, &
+ phi)
+ ! add metadata for reassembly
+ flavor1 = myrank
+ flavor2 = iglob
+ write(IOUT,*) theta, phi, radius, flavor1, flavor2
+ enddo
+
+ enddo
+
+ close(IOUT)
+
+end program write_CitcomS_tracers
More information about the cig-commits
mailing list