Ignore:
Timestamp:
Aug 5, 2008, 2:26:34 PM (16 years ago)
Author:
ole
Message:

Fixed urs2sts such that start and stop times are updated across mux sources

Location:
anuga_core/source/anuga/shallow_water
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5606 r5611  
    49034903
    49044904    # Call underlying C implementation urs2sts_ext.c   
    4905     data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
     4905    data = read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
    49064906
    49074907    msg='File parameter values were not read in correctly from c file'
     
    49454945        quantity[i][:]=data[i][:parameters_index] # Was data[i][:-OFFSET]
    49464946       
     4947        first_time_step = data[i][parameters_index+3]
     4948        #print 'datamanager:', i, first_time_step, dt*first_time_step
    49474949        starttime=min(dt*data[i][parameters_index+3],starttime)
    49484950       
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5602 r5611  
    65716571        time_start_e = time_start_z
    65726572
     6573        # FIXME(Ole): Should they really be treated differently?
    65736574        time_start_n = array([[10.0,11.5,13,14.5,17.7],
    65746575                              [10.6,11.8,12.9,14.2,17.4],
     
    68126813        msg = 'sts starttime was %f. Should have been %f'\
    68136814            %(sts_starttime, starttime)
    6814         assert allclose(sts_starttime-delta_t, starttime), msg
    6815         #assert allclose(sts_starttime, starttime), msg
     6815        #assert allclose(sts_starttime-delta_t, starttime), msg
     6816        assert allclose(sts_starttime, starttime), msg
    68166817        ## FIXME - have done a dodgy to get it through here ###   
    68176818   
     
    1019910200    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
    1020010201    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_individual_sources')
    10201     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_combined_sources')     
     10202    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts')       
    1020210203
    1020310204   
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5602 r5611  
    5757        }
    5858    }
     59   
    5960    if (nft[ista] != -1)
    6061    {
     
    6869        }
    6970    }     
     71   
    7072    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
    7173    {
    7274        /* gauge never started recording, or was outside of all grids,
    73     fill array with 0 */
     75        fill array with 0 */
    7476        for(it = 0; it < nt; it++)
    7577        {
     
    9193            /* deal with the tide gauge at hand */
    9294            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
     95            {
    9396                /* gauge is recording at this time */
    94             {
    9597                memcpy(data + it, muxData + offset, sizeof(float));
    9698                offset++;
     
    116118        if(last_it < nt - 1)
    117119            /* the loop was exited early because the gauge had
    118         finished recording */
     120            finished recording */
    119121            for(it = last_it+1; it < nt; it++)
    120122                data[it] = NODATA;
     
    267269
    268270        /* Read the start and stop times for this source */
    269         fread(fros + i*(*total_number_of_stations), *total_number_of_stations*sizeof(int), 1, fp);
    270         fread(lros + i*(*total_number_of_stations), *total_number_of_stations*sizeof(int), 1, fp);
     271        fread(fros + i*(*total_number_of_stations),
     272              *total_number_of_stations*sizeof(int), 1, fp);
     273        fread(lros + i*(*total_number_of_stations),
     274              *total_number_of_stations*sizeof(int), 1, fp);
    271275
    272276        /* Compute the size of the data block for this source */
     
    291295
    292296    // Store time resolution and number of timesteps   
    293     // These are the same for all stations, so
     297    // These are the same for all stations as tested above, so
    294298    // we take the first one.
    295299    *delta_t = (double)mytgs0[0].dt;
     
    312316    FILE *fp;
    313317    int total_number_of_stations, i, isrc, ista, k;
    314     //struct tgsrwg* mytgs0=NULL;
    315318    char *muxFileName;
    316     int istart, istop;
     319    int istart=-1, istop=-1;
    317320    int number_of_selected_stations;
    318321    float *muxData=NULL; // Suppress warning
     
    325328    long int offset;
    326329
    327     int number_of_time_steps;
     330    int number_of_time_steps, N;
    328331    double delta_t;
    329     //long numDataMax;
     332   
     333    // Shorthands pointing to memory blocks for each source
     334    int *fros_per_source=NULL;     
     335    int *lros_per_source=NULL;         
    330336   
    331337    _read_mux2_headers(numSrc,
     
    370376
    371377    // For each selected station, allocate space for its data
    372 
    373     len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
     378    len_sts_data = number_of_time_steps + POFFSET; // Max length of each timeseries?
    374379    for (i = 0; i < number_of_selected_stations; i++)
    375380    {
     
    381386            return NULL;
    382387        }
    383 
    384         ista = (int) permutation[i]; // Get global index into mux data
    385    
    386         sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
    387         sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
    388         sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
    389         sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
    390         sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
    391388    }
    392389
     
    394391
    395392    muxData = (float*)malloc(numDataMax*sizeof(float));
     393   
    396394    /* Loop over all sources */
    397     //FIXME: remove istart and istop they are not used.
    398     istart = -1;
    399     istop = -1;
    400395    for (isrc = 0; isrc < numSrc; isrc++)
    401396    {
     397   
     398        // Shorthands to local memory
     399        fros_per_source = (int*) fros + isrc*total_number_of_stations;
     400        lros_per_source = (int*) lros + isrc*total_number_of_stations;     
     401           
     402   
     403   
    402404        /* Read in data block from mux2 file */
    403405        muxFileName = muxFileNameArray[isrc];
     
    415417        fseek(fp, offset, 0);
    416418
    417         numData = getNumData(fros + isrc*total_number_of_stations,
    418                              lros + isrc*total_number_of_stations,
     419        numData = getNumData(fros_per_source,
     420                             lros_per_source,
    419421                             total_number_of_stations);
    420422                             
     
    431433            ista = (int) permutation[i]; // Get global index into mux data 
    432434       
    433             /* fill the data0 array from the mux file, and weight it */
     435            // fill the data0 array from the mux file, and weight it
    434436            fillDataArray(ista,
    435437                          total_number_of_stations,
    436                           mytgs0[ista].nt,
     438                          number_of_time_steps,
    437439                          mytgs0[ista].ig,
    438                           fros + isrc*total_number_of_stations,
    439                           lros + isrc*total_number_of_stations,
     440                          fros_per_source,
     441                          lros_per_source,
    440442                          temp_sts_data,
    441443                          &istart,
     
    443445                          muxData);
    444446
    445             /* weight appropriately and add */
     447            // Weight appropriately and add
    446448            for(k = 0; k < mytgs0[ista].nt; k++)
    447449            {
     
    455457                }
    456458            }
     459           
     460            // Update metadata (e.g. start time and end time)
     461            N = number_of_time_steps;
     462           
     463           
     464            /*
     465            printf("station =  %d, source = %d, fros = %f, lros = %f\n",
     466                   ista,
     467                   isrc,
     468                   (float) fros_per_source[ista],
     469                   (float) lros_per_source[ista]);                 
     470            */
     471                   
     472                   
     473                   
     474            if (isrc == 0) {
     475                // Assign values for first source
     476                sts_data[i][N] = (float)mytgs0[ista].geolat;
     477                sts_data[i][N+1] = (float)mytgs0[ista].geolon;
     478                sts_data[i][N+2] = (float)mytgs0[ista].z;
     479                sts_data[i][N+3] = (float)fros_per_source[ista];
     480                sts_data[i][N+4] = (float)lros_per_source[ista];
     481            } else {
     482                // Update first and last timesteps for subsequent sources
     483                if (sts_data[i][N+3] > (float)fros_per_source[ista]) {         
     484                    if (verbose) {
     485                        printf("Adjusting start time for station %d and source %d",
     486                               ista, isrc);
     487                        printf(" from %f to %f\n", sts_data[i][N+3], (float)fros_per_source[ista]); 
     488                    }
     489                    sts_data[i][N+3] = (float)fros_per_source[ista];
     490                }
     491               
     492               
     493                if (sts_data[i][N+4] < (float)lros_per_source[ista]) {         
     494                    if (verbose) {
     495                        printf("Adjusting end time for station %d and source %d",
     496                               ista, isrc);
     497                        printf(" from %f to %f\n", sts_data[i][N+4], (float)lros_per_source[ista]); 
     498                    }
     499                    sts_data[i][N+4] = (float)lros_per_source[ista];
     500                }               
     501            }
    457502        }
    458503    }
Note: See TracChangeset for help on using the changeset viewer.