main.c

用于2维的射线追踪

/*
 * Ray2mesh : software for geophysicists.
 * Compute various scores attached to the mesh cells, based on geometric
   information that rays bring when the traverse the cell.
 *
 * Copyright (C) 2003, St閜hane Genaud and Marc Grunberg
 *
 * This tool is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free
 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#ifdef USE_MPI
#include <mpi.h>
#endif
#ifdef USE_ZLIB
#include "zlib.h"
#endif
#include <const_ray2mesh.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <sys/time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <errno.h>
#include <signal.h>

#include <popt.h>


#include <iasp/iasp.h>
#include <ray/raydescartes.h>

#include <mesh/mesh.h>
#include <mesh/layer.h>
#include <mesh/cellinfo.h>
#include <mesh/mesh2xml.h>

#include <sparse/sparse.h>

#include "ray2mesh.h"
#include "r2mfile.h"
#include "merge.h"
#include "score.h"
#include "buffer.h"
#include "filter.h"


char *celldatafilename = NULL;		/**< these are kept global to close */
FILE *fdinput = NULL;			/**< them in case of emergency (Ctrl-C) */

struct mesh_t *mesh;			/**< mesh structure */
struct cell_info_t ****cell_info = NULL;/**< traversed cells information */

char *output_format = NULL;			/**< r2m and/or sco output */
int rank;

volatile int mem_used=0;

/** \brief shows ray2mesh version ans compilation flags */
void ray2mesh_info(struct ray_config_t *ray_config)
{
    char *modelver = libmodelversion();
    char *rayver = librayversion();
    char *raycodever = libraycodeversion();
    char *meshver = libmeshversion();
    char *sparsever = libsparseversion();

    fprintf(stderr, "%s: version %s", PACKAGE, VERSION);
#ifdef STATIC_BIN
    fprintf(stderr, " [*static binary*]");
#endif
    fprintf(stderr, "\nusing \t%s\n", modelver);
    fprintf(stderr, "\t%s\n", rayver);
    fprintf(stderr, "\t%s\n", raycodever);
    fprintf(stderr, "\t%s\n", meshver);
    fprintf(stderr, "\t%s\n", sparsever);

    free(modelver);
    free(rayver);
    free(raycodever);
    free(meshver);
    free(sparsever);

    fprintf(stderr, "Compilation flags : ");

#ifdef USE_MPI
    fprintf(stderr, "USE_MPI ");
    #ifdef MASTER_SLAVE
        fprintf(stderr, "MASTER_SLAVE ");
    #else
        fprintf(stderr, "SCATTER ");
    #endif

    #ifdef USE_MPILB
    fprintf(stderr, "USE_MPILB ");
    #endif
#endif

#ifdef DEBUG
    fprintf(stderr, "DEBUG ");
#endif
#ifdef ALLOC_MEMORY_CHUNK
    fprintf(stderr, "ALLOC_MEMORY_CHUNK ");
#endif
#ifdef RAYTRACING_ONLY
    fprintf(stderr, "RAYTRACING_ONLY ");
#endif
#ifdef USE_ZLIB
    fprintf(stderr, "USE_ZLIB(%s) ", zlibVersion());
#endif
    fprintf(stderr, "\n");
}

/** \brief given a filename, stops if file is unreadable  */
void check_file_access(char *filename)
{
struct stat st;
    /* check read permissions */
    if (access(filename, R_OK) < 0) {
	fprintf(stderr, "%s : cannot access file %s. Exiting.\n", PACKAGE, filename);
	exit(1);
    }
    /* check filetype */
    stat(filename,&st);
    if (!S_ISREG(st.st_mode)) {
	fprintf(stderr, "%s : file %s is not a regular file. Exiting.\n", PACKAGE, filename);
	exit(1);
    }
}

/** \brief handler for Ctrl-C, close the files */
void emergency_halt()
{
    fprintf(stdout, "*** Emergency stop (Ctrl-C). Sync files ***\n");
    fflush(stdout);

    if (strchr(output_format, 'r')) {
	/* R2M */
	make_domain_info_file("emergency.r2m", cell_info, mesh, 0, 1, 0);
    } else {
	/* SCO */
        fprintf(stdout, "-> Computing cell score ... ");
	compute_score(cell_info, mesh);
        fprintf(stdout, "done.\n-> Dumping cell information ... ");
	mesh_cellinfo_write_sco ("emergency.sco", cell_info, mesh);
    }
    fprintf(stdout, "done.\n");
#ifdef USE_MPI
    MPI_Finalize();
#endif
    exit(1);
}

/*
 * \brief parse command line using popt : checks validity of optoins 
 * specified and returns parameters values.
 *
 * PRECOND : the default values for optional arguments must be set.     
 */
void
parse_command_line(int argc, char **argv,
		   struct ray_config_t *ray_config,
		   struct ray_filter_t *ray_filter,
		   char **meshfile,
		   char **outputfile,
		   char **inputdatafile,
		   char **rayfilter_filename, 
		   float *mem_limit, 
		   char **output_format,
		   char **tmpdir
#ifdef MASTER_SLAVE
		   , int *chunksize
#endif
    )
{
    static float length_opt;
    static char *meshfile_opt=NULL;	/* there to build on IRIX with popt */
    static char *outputfile_opt=NULL;
    static char *vmodel;
    static char *inputdatafile_opt=NULL;
    static char *iasptables_opt;
    static char *filtered=NULL;
    static float mem_limit_opt;
    static char *setfilter_opt=NULL;
    static char *outputformat_opt=NULL;
    static int iterate_opt;
    static char *tmpdir_opt=NULL;

#ifdef MASTER_SLAVE
    static int chunksize_opt;
#endif
    
    poptContext cx;
    int rc;
    int nbarg;
    

#include "options.h"

    /* init */
    set_ray_config_to_default(ray_config);

    /* ray_config->iterative_mode = ITERATIVE_MODE_ANGULAR; */
    /* ray_config->iterative_mode = ITERATIVE_MODE_OFF; */
    /* ray_config->iterative_mode = ITERATIVE_MODE_ON; */

    /* args parsing using popt */
    cx = poptGetContext(PACKAGE, argc, (const char **) argv, options, 0);
    if (argc == 1) {
	poptPrintUsage(cx, stdout, 0);
	exit(1);
    }

    do {
	rc = poptGetNextOpt(cx);
	switch (rc) {
	    /* help */
	case OPT_HELP:
	    poptPrintHelp(cx, stdout, 0);
	    exit(0);
	case OPT_VELOCITY_MODEL:
            check_file_access(vmodel);
	    ray_config->velocity_model = load_velocity_model(vmodel);
	    ray_config->get_velocity_p = my_get_velocity_p;
	    ray_config->get_velocity_s = my_get_velocity_s;
	    break;
	case OPT_LENGTH:
	    ray_config->length_inc = length_opt;	/* update if increment specified */
	    break;
	case OPT_OUTPUT_FORMAT:	/* sco output selected */
	    if ( !strchr(outputformat_opt,'s') &&
                 !strchr(outputformat_opt,'r'))
            {
                    fprintf(stderr, "format must be 'sco' or 'r2m'.\n");
                    exit(1);
            }
	    if (strchr(outputformat_opt,'s') && strchr(outputformat_opt,'r')) {
	            fprintf(stderr, "sorry, both output not yet implemented\n");
                    exit(1);
	    }
	    break;
#ifdef MASTER_SLAVE
	case OPT_CHUNKSIZE:
	    *chunksize = chunksize_opt;
	    break;
#endif
	case OPT_TABLES:
#ifdef DEBUG
	    fprintf(stderr, "setting iasp tables path to \"%s\"\n",
		    iasptables_opt);
#endif
	    set_modnam(iasptables_opt);
	    break;
	case OPT_ITERATE:
	    ray_config->iterative_mode = ITERATIVE_MODE_ON;
	    break;
	case OPT_MEM_LIMIT:
	    *mem_limit = mem_limit_opt;
	    break;
	case OPT_SET_FILTER:
	    nbarg = parse_filter_value(strdup(setfilter_opt), ray_filter);
	    ray_filter->activated = 1;
	    if (nbarg!=3) {
		    fprintf(stderr, "parsing filter string (%s) failed\n", 
				    setfilter_opt); 
		    fprintf(stderr, "use --set-filter=r,d,n\n");
		    exit(1);
	    }
	    break;
	case OPT_VERSION:
	    ray2mesh_info(ray_config);
	    exit(0);
	case POPT_ERROR_BADOPT:	/* error */
	    fprintf(stderr, "%s: %s\n", poptBadOption(cx, 0),
		    poptStrerror(rc));
	    poptPrintUsage(cx, stdout, 0);
	    exit(1);
	}
    }
    while (rc > 0);

    /* check meshfile */
    if (!meshfile_opt) {
	fprintf(stderr, "%s: missing -m option\n", PACKAGE);
	exit(1);
    } else {
	*meshfile = strdup (meshfile_opt);
    }

    if (!outputfile_opt) {
        fprintf(stderr, "%s: missing -o option\n", PACKAGE);
	exit(1);
    } else {
    	*outputfile = strdup(outputfile_opt);
    }

    /* check inputdatafile option */
    if (!inputdatafile_opt) {
	fprintf(stderr, "%s: missing -i option\n", PACKAGE);
	exit(1);
    } else {
    	*inputdatafile = strdup(inputdatafile_opt);
    }

    /* check output format */
    if (!outputformat_opt) {
        fprintf(stderr, "%s: missing -f option, sco or r2m\n", PACKAGE);
	exit(1);
    } else {
    	*output_format = strdup(outputformat_opt);
    }

    if (filtered) {
    	*rayfilter_filename = strdup (filtered);
    }

    /* set tmp directory */
    if (tmpdir_opt) {
    	*tmpdir = strdup (tmpdir_opt);
    } else {
	*tmpdir = strdup ("/tmp");
    }
	
    {
     struct stat stabuf;

     if (stat(*tmpdir, &stabuf) == -1) {
        perror(*tmpdir);
	exit(1);
     }
     if (!S_ISDIR(stabuf.st_mode)) {
	fprintf(stderr, "%s is not a directory !\n", *tmpdir);
	exit(1);
     }
    }

    poptFreeContext(cx);
}



void print_ray_ko_info (int rank, struct raydata_t *raydata, int i, int nbrays, char *rai_status) {
	fprintf(stderr,
		   "[process %d] KO %d/%d,  %s,  -s %f,%f,%f -d %f,%f,%f -p %s\n",
		   rank,
		   i,nbrays,
		   rai_status?rai_status:"",
		   raydata[i].src.lat * TO_DEG,
		   raydata[i].src.lon * TO_DEG,
		   raydata[i].src.prof,
		   raydata[i].dest.lat * TO_DEG,
		   raydata[i].dest.lon * TO_DEG,
		   raydata[i].dest.prof, raydata[i].phase); 
}


/* \brief given a set of raydata elements (src and dest) pairs 
 * for rays, and the number of rays contained in the set,               
 * run ray3D_descartes of this bunch of ray data with a discretization  
 * step length of length_inc km. nb_ray_computed and nb_ray_rejected are
 * updated accordingly                                                  
 * 
 * POSTCOND : raydata is freed before returning                         
 */
void bunch_of_ray(struct raydata_t *raydata,
		  const int nbrays,
		  const long int offset,
		  struct ray_config_t *ray_config,
		  struct ray_filter_t *ray_filter,
		  int *nb_ray_computed,
		  int *nb_ray_rejected,
		  FILE * filterfd,
		  FILE * matrix_length_fd, 
		  FILE * residual_time_fd,
		  FILE * event_fd)
{
    struct ray3D_t **my3Drays, **ptr;
    int i;
    char **phases, **cp_phases;
    char *rai_status=NULL;
    double dT;			/* time residual */
    double delta;

    fprintf(stdout, "*** [process %d] bunch_of_ray rays  rayId=[%d  --> %d]\n",
	    rank, *nb_ray_computed, *nb_ray_computed + nbrays -1);



    for (i = 0; i < nbrays; i++) {

	if (ray_filter->activated) {
	        /* check if nb bundle is ok */
		if (raydata[i].event_id < ray_filter->nb_event_min) {
			print_ray_ko_info (rank, raydata, i, nbrays, "nbbundle");
			(*nb_ray_rejected)++;
			continue;
		}
	        /* check if delta is ok */
		delta = AOB(&raydata[i].src, &raydata[i].dest)*TO_DEG;
		if (delta > ray_filter->delta_max) {
			print_ray_ko_info (rank, raydata, i, nbrays, "delta");
			(*nb_ray_rejected)++;
			continue;
		}
	}

	/* do the ray-tracing */
	phases = cp_phases = parse_phase(raydata[i].phase);
	my3Drays = ray3D_descartes(phases, &raydata[i].src, 
			&raydata[i].dest, ray_config);

	/* desalloc phases */
	while (*cp_phases) {
	    free(*cp_phases);
	    cp_phases++;
	}
	free(phases);

	/* loop over all returned rays phases for this ray */
	for (ptr = my3Drays; *ptr; ptr++) {
	    /* dT = Tobserved - Tcomputed */
	    dT = raydata[i].ray_travel_time - (*ptr)->time;

	    /* check if residual time is small enough */
	    if (ray_filter->activated && 
	        fabs(dT) > ray_filter->residual_max) 
	    {
	    	/* free */
		print_ray_ko_info (rank, raydata, i, nbrays, "residual");
            	(*nb_ray_rejected)++;
            	free_ray3D_descartes(*ptr);
            	*ptr = NULL;
		continue;
	    }
	
	    rai_status = str_ray_status((*ptr)->status);

	    /* ray-trace failed */
	    if ((*ptr)->status != RAYTRACE_OK) {
		print_ray_ko_info (rank, raydata, i, nbrays, rai_status);
		(*nb_ray_rejected)++;
		free_ray3D_descartes(*ptr);
            	*ptr = NULL;
		free(rai_status);
		rai_status=NULL;
		continue;
	    }

	    /* ray-trace is ok */
	    fprintf(stderr, "[process %d] OK ray %d/%d (phase %s) %s length=%f nbiter=%d dist_err=%f\n",
		    rank,
		    i, nbrays, raydata[i].phase, rai_status, 
		    (*ptr)->total_length,
		    (*ptr)->nb_iter_angle+(*ptr)->nb_iter_length,
		    (*ptr)->dist_err);
	    free(rai_status);
	    rai_status=NULL;

	    /* ray is ok */
	    if (filterfd) {
		/* user want to keep info about good rays */
		fprintf(filterfd, "%ld %g %g %g %g %g %g %s %g %f\n",
			raydata[i].event_id,
			raydata[i].src.lat * TO_DEG,
			raydata[i].src.lon * TO_DEG,
			raydata[i].src.prof,
			raydata[i].dest.lat * TO_DEG,
			raydata[i].dest.lon * TO_DEG, 
			raydata[i].dest.prof,
			raydata[i].phase,
			/*(*ptr)->code, */
			raydata[i].ray_travel_time, 
			dT);
	    }

#ifndef RAYTRACING_ONLY
	    /* we must compute the time residual */
	    if (residual_time_fd) {
		fprintf(residual_time_fd, "%ld %lf\n",
			(long int) *nb_ray_computed, (double) dT);
	    }

	    /* write event_id ray_id travel_time */
   	    if (event_fd) {
		fprintf(event_fd, "%ld %ld %lf\n", 
			raydata[i].event_id, (long int) *nb_ray_computed, 
			(*ptr)->time);
	    }

	    /* pass previous cell_info to accumulate results 
	     * from several rays */
	    cell_info = ray2mesh2(*ptr, *nb_ray_computed, mesh, cell_info,
			 matrix_length_fd);
#endif

	    /* free */
	    free_ray3D_descartes(*ptr);
	    *ptr = NULL;
	    (*nb_ray_computed)++;
	}

	free(my3Drays);
    }

    free(raydata);
    
}


/*-----------------------------------------*/
/* main                                    */
/*------- ---------------------------------*/
#ifdef USE_MPI
#ifdef MASTER_SLAVE
#include "main_mpi_ms.c"
#else
#include "main_mpi.c"
#endif
#else
#ifdef DEBUG
#include <mcheck.h>
#endif
#include "main_seq.c"
#endif