main_seq.c

用于2维的射线追踪

/*
 * Ray2mesh : software for geophysicists.
 * Compute various scores attached to the mesh cells, based on geometric
   information that rays bring when they traverse cells.
 *
 * 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.
 */

#include<stdlib.h>
#include<string.h>

#include "save_to_disk.h"
#include "memusage.h"
#include "renumber.h"

#define  DEFAULT_SLICE 50			/**< the number of raydata to distribute in one round */

/*----------------------------------------------------------------------*/
/* main (sequential code)                                               */
/*----------------------------------------------------------------------*/
int main(int argc, char **argv)
{				/* return the nb of ray computed */

    struct raydata_t *raydata;
    struct ray_config_t ray_config;
    struct ray_filter_t ray_filter;
    char *meshfile = NULL;
    char *inputdatafile = NULL;
    char *out_filename = NULL;		/**< cell_data filename really used (may change with --limit */
    char *tmpdir=NULL;			/**< only usefull with mpi version */
    char *rayfilter_filename = NULL;
    FILE *filter_fd = NULL;
    
    int nb_ray_total = 0;
    int nb_ray_computed = 0;
    int nb_ray_rejected = 0;
    
    int iterate = 0;			/**< is it a unique ray (0) or an input datafile (1) */
    int nbread = 0;
    int nberr = 0;
    
    float mem_used_start=0;		/**< memory used at the begin of the raytracing process */

    float limit = -1.;
    float memsize = 0.;
    int slice = 0;			/**< size of ray to read */
    int nb_ray_stored = 0;		/**< nb of ray used to dump memory to disk */
    int nb_save = 0;	   		/**< nb of times we have dumped memory to disk */

    int size;				/**< nb of ray to trace */

    char *xml_output;
    
    FILE *sparse_fd;
    FILE *res_fd;
    FILE *event_fd;

    /* catch Ctrl-C signal */
    signal(SIGINT, emergency_halt);

#ifdef DEBUG
    /* remember to set MALLOC_TRACE=file */
    mtrace();
#endif

    /******************/
    /* parse cmd line */
    /******************/
    ray_filter.activated=0;
    parse_command_line(argc, argv,
		       &ray_config,
		       &ray_filter,
		       &meshfile, 
		       &out_filename, 
		       &inputdatafile, 
		       &rayfilter_filename, 
		       &limit, 
		       &output_format,
		       &tmpdir);

    check_file_access(meshfile);
    check_file_access(inputdatafile);

    fprintf(stderr, "Raytracing mode is ");   
    switch (ray_config.iterative_mode) {
	case ITERATIVE_MODE_OFF:
	    fprintf(stderr, "ITERATIVE_MODE_OFF\n");
	    break;
	case ITERATIVE_MODE_ANGULAR:
	    fprintf(stderr, "ITERATIVE_MODE_ANGULAR\n");
	    break;
	case ITERATIVE_MODE_ON:
	    fprintf(stderr, "ITERATIVE_MODE_ON\n");
    }

    if (ray_filter.activated) { 
     fprintf(stderr, "Filter activated, removing rays with :\n");
     fprintf(stderr, "\tresidual  > %.2fs\n\tdist epi  > %.2fdeg\n\tnb bundle > %d\n\n",
		     ray_filter.residual_max, 
		     ray_filter.delta_max, 
		     ray_filter.nb_event_min);
    } else {
	fprintf(stderr, "Filter disable\n");
    }


    /* open the filter file if needed */
    if (rayfilter_filename) {
       char *filename;
       filename = (char *) malloc(sizeof(char)*(
                   strlen(rayfilter_filename) + strlen(".sei") + 1));
       assert(filename);
       sprintf(filename, "%s.sei", rayfilter_filename);

       if ((filter_fd = fopen(filename, "w")) == NULL) {
          fprintf(stderr, "Cannot open %s for writing ray filtered.\n", 
                          filename);
	  exit(1);
       }
       free(filename);
    } 



    /***********************************************/
    /* set chunk size value ie. nb of rays to send */
    /***********************************************/
    if (limit > 0) {
       /* if user specified a memory limit */
       /* slice size if 1/10 of the *estimated* value */
       mem_used_start = get_mem_usage();
       slice = limit * 1024 / (10 * ONE_RAY_MEM);

       fprintf(stdout, "*** limit memory usage to %.1fM\n", limit);
       fprintf(stdout, "*** memory limit specified ... forcing slice to %d rays.\n",
			slice);
    } else {
	slice = DEFAULT_SLICE;
	fprintf(stdout, "*** setting slice to %d rays.\n", slice);
    }
    fflush(stdout);

    

    /* init of the mesh->parameters struct  */
    if (!(mesh = mesh_init_from_file(meshfile))) {
	exit(1);
    }

    /* open rays input file */
    if (!(fdinput = fopen(inputdatafile, "r"))) {
	fprintf(stderr, "Can not open input ray file '%s' ... exiting!\n",
			inputdatafile);
	exit(1);
    }
	
    fprintf(stdout, "** Using data file '%s' for ray computing **\n\n", 
		    inputdatafile);
    size = get_number_of_lines(fdinput);
    


    /***************************************/
    /* open the sparse, res and event file */
    /***************************************/
    sparse_fd = NULL;
    res_fd = NULL;
    event_fd = NULL;
    change_to_next_files (mesh, out_filename, 
		    	  size, 
		          &sparse_fd, &res_fd, &event_fd, 
		          0 /*rank*/ , 
			  0 /*nb_save*/, 
			  1 /*want_sparse_file*/);

    /*************************************/
    /* iterate over input datafile lines */
    /*************************************/
    while (1) {
	/* read a bunch of lines from input file */
	raydata = get_raydata(fdinput, slice, &nbread, &nberr);
	nb_ray_total += (nbread + nberr);



	fprintf(stdout, "\n*** requesting SLICE=%d / nbread=%d\n",
		slice, nbread);

	if (!raydata) {
	    break;
	}

	bunch_of_ray(raydata, nbread, 0,	/* Put 0 as offset */
		     &ray_config, &ray_filter,
		     &nb_ray_computed, &nb_ray_rejected,
		     filter_fd, sparse_fd, res_fd, event_fd);

	fprintf(stdout,
		"*** computed %d/%d rays and rejected %d/%d (%.2f%%)\n",
		nb_ray_computed, nb_ray_total, nb_ray_rejected,
		nb_ray_total,
		(float) nb_ray_rejected / (float) (nb_ray_total) * 100);

	/* shows memory usage */
	if (limit > 0) {
	    memsize = get_mem_usage() - mem_used_start;
	    fprintf(stdout, "*** MEM(pid=%d) = %.2fMB\n", getpid(), memsize);
	}

	iterate = !feof(fdinput);
	if (!iterate) {
	    fclose(fdinput);
	    break;
	}

	/* save memory to disk */
	if (    (limit > 0) 
             && (memsize > limit)
	     && (nb_ray_computed != nb_ray_stored))  
	{
             save_memory_to_disk(output_format, out_filename, cell_info, mesh, 
			     0 /*rank*/, 1 /*nbprocs*/, nb_save);
	     nb_ray_stored = nb_ray_computed;
	     nb_save++;     

	     /* next set of files */
	     change_to_next_files (mesh, out_filename, size, 
			           &sparse_fd, &res_fd, &event_fd, 
				   0 /* rank */, 
				   nb_save, 
				   1 /*want_sparse_file*/);

	     fprintf(stderr, "\n*** memory limit (%f MB) reached after %d rays.\n", 
			     limit, nb_ray_computed);
	}
    }


    /*********************************************/
    /* save the data remaining in memory to disk */
    /*********************************************/
    if (    (limit > 0) 
	 && (nb_ray_computed != nb_ray_stored)) 
    { 
	    save_memory_to_disk(output_format, out_filename, cell_info, mesh, 
			         0 /*rank*/, 1 /*nbprocs*/, nb_save); 
    } 
    
    /* close sparse, res and evt files */
    fclose(sparse_fd);
    fclose(res_fd);
    fclose(event_fd);

    /* everyone close its filtered ray file (if filtering used) */
    if (filter_fd) {
	fclose(filter_fd);
    }
	
    fprintf(stdout,
	    "*** computed %d/%d rays and rejected %d/%d (%.2f%%)****\n",
	    nb_ray_computed, nb_ray_total, nb_ray_rejected, nb_ray_total,
	    (float) nb_ray_rejected / (float) (nb_ray_total) * 100);

#ifndef RAYTRACING_ONLY
    /**************************/
    /* output r2m or sco file */
    /* all data is in memory  */
    /**************************/
    if (limit < 0)  {
       if (strchr(output_format,'r')) {
	    char *r2m_file;
	    
	    /* r2m */
	    r2m_file = construct_filename(out_filename, "r2m", 0, 0);
	    fprintf(stdout,
		    "*** writing cell data (r2m formated) in %s\n", r2m_file);
	    mesh_add_data_filename(mesh, R2M, r2m_file);
	    make_domain_info_file(r2m_file, cell_info, mesh, 0, 1, 0);
	    free(r2m_file);
       } else {
	    /* sco */
	    char *sco_file;

    	    fprintf(stdout, "*** score computing ***\n");
    	    compute_score(cell_info, mesh);

	    sco_file = construct_filename(out_filename, "sco", 0, 0);
	    fprintf(stdout,
		    "*** writing cell data (sco formated) in %s\n", sco_file);
	    mesh_add_data_filename(mesh, SCO, sco_file);
	    mesh_cellinfo_write_sco (sco_file, cell_info, mesh);
	    free(sco_file);
       }
    }


    /* renumber sparse file only if all data remains in memory */
    /* if not you should use the merger software               */
    if (limit < 0) {
       char * filename;
       filename = mesh->data[SPARSE]->filename[0];
       renumber_sparse_file (0, nb_ray_computed, filename);

       filename = mesh->data[RES]->filename[0];
       renumber_res_file (0, nb_ray_computed, filename);

       filename = mesh->data[EVT]->filename[0];
       renumber_evt_file (0, nb_ray_computed, filename);
    }
    

    /* save the xml enrichied with SCO/R2M/SPARSE/RES sections */
    xml_output = (char *) malloc((strlen(out_filename) + strlen(".xml") + 1) *
			sizeof(char));
    assert(xml_output);
    sprintf(xml_output, "%s.xml", out_filename);
    mesh2xml(mesh, xml_output);
    free(xml_output);

#endif
    
    free_velocity_model (ray_config.velocity_model);


#ifdef DEBUG
    muntrace();
#endif

    return (0);
}