Commit 18ceef00 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

samplegridicon bug fix.

parent fb0915df
......@@ -7,7 +7,6 @@ typedef struct {
int ncells;
int *neighbor; // neighbor cell index
int *parent; // parent cell index
int *data; // field data index for parent
int *child; // child cell index
const char *filename;
} cellindex_type;
......@@ -25,7 +24,6 @@ void free_cellindex(cellindex_type *cellindex)
{
if ( cellindex->neighbor ) Free(cellindex->neighbor);
if ( cellindex->parent ) Free(cellindex->parent);
if ( cellindex->data ) Free(cellindex->data);
if ( cellindex->child ) Free(cellindex->child);
Free(cellindex);
}
......@@ -72,10 +70,9 @@ cellindex_type *read_cellindex(const char *filename)
cellindex_type *cellindex = (cellindex_type*) Malloc(sizeof(cellindex_type));
cellindex->ncells = ncells;
cellindex->neighbor = (int*) Malloc(3*ncells*sizeof(int));
cellindex->neighbor = NULL;
// cellindex->neighbor = (int*) Malloc(3*ncells*sizeof(int));
cellindex->parent = (int*) Malloc( ncells*sizeof(int));
cellindex->data = (int*) Malloc( ncells*sizeof(int));
for ( int i = 0; i < ncells; ++i ) cellindex->data[i] = i;
cellindex->child = NULL;
// cellindex->child = (cid != CDI_UNDEFID) ? (int*) Malloc(4*ncells*sizeof(int)) : NULL;
double *data = (double *) Malloc(ncells*sizeof(double));
......@@ -85,18 +82,18 @@ cellindex_type *read_cellindex(const char *filename)
{
int varID, levelID, nmiss;
streamInqRecord(streamID, &varID, &levelID);
if ( varID == nid || varID == pid /*|| varID == cid*/ )
if ( varID == pid /* || varID == nid || varID == cid */ )
{
streamReadRecord(streamID, data, &nmiss);
if ( varID == nid ) copy_data_to_index(ncells, data, cellindex->neighbor+levelID*ncells);
else if ( varID == pid ) copy_data_to_index(ncells, data, cellindex->parent);
// else if ( varID == cid ) copy_data_to_index(ncells, data, cellindex->child+levelID*ncells);
if ( varID == pid ) copy_data_to_index(ncells, data, cellindex->parent);
// else if ( varID == nid ) copy_data_to_index(ncells, data, cellindex->neighbor+levelID*ncells);
// else if ( varID == cid ) copy_data_to_index(ncells, data, cellindex->child+levelID*ncells);
}
}
// Fortran to C index
for ( int i = 0; i < 3*ncells; ++i ) cellindex->neighbor[i] -= 1;
for ( int i = 0; i < ncells; ++i ) cellindex->parent[i] -= 1;
// for ( int i = 0; i < 3*ncells; ++i ) cellindex->neighbor[i] -= 1;
streamClose(streamID);
......@@ -172,7 +169,7 @@ int find_index(int search, int n, const int *restrict array)
typedef struct
{
int p, d;
int p, i;
}
sinfo_t;
......@@ -194,21 +191,23 @@ void compute_child(cellindex_type *cellindex1, cellindex_type *cellindex2)
{
int ncells1 = cellindex1->ncells;
int *parent1 = cellindex1->parent;
int *idx1 = (int*) Malloc(ncells1*sizeof(int));
for ( int i = 0; i < ncells1; ++i ) idx1[i] = i;
for ( int i = 1; i < ncells1; ++i )
if ( parent1[i] < parent1[i-1] )
{
if ( cdoVerbose ) cdoPrint("Parent index not sorted in %s (index %d)!", cellindex1->filename, i);
if ( cdoVerbose ) cdoPrint("Sort parent index of %s!", cellindex1->filename);
sinfo_t *sinfo = (sinfo_t*)Malloc(ncells1*sizeof(sinfo_t));
for ( int j = 0; j < ncells1; ++j )
{
sinfo[j].p = parent1[j];
sinfo[j].d = cellindex1->data[j];
sinfo[j].i = idx1[j];
}
qsort(sinfo, ncells1, sizeof(sinfo_t), cmpsinfo);
for ( int j = 0; j < ncells1; ++j )
{
parent1[j] = sinfo[j].p;
cellindex1->data[j] = sinfo[j].d;
idx1[j] = sinfo[j].i;
}
Free(sinfo);
break;
......@@ -225,10 +224,12 @@ void compute_child(cellindex_type *cellindex1, cellindex_type *cellindex2)
for ( int k = 0; k < 4; ++k )
{
if ( i != parent1[j+k] ) break;
child2[i*4+k] = j+k;
// child2[i*4+k] = j+k;
child2[i*4+k] = idx1[j+k];
}
// if ( i%10000 == 0 ) printf("%d %d %d %d %d %d\n", i, j, parent1[j], parent1[j+1], parent1[j+2], parent1[j+3]);
}
Free(idx1);
}
static
......@@ -244,7 +245,6 @@ void compute_sum(int i, int *n, double *sum, double *sumq, int kci, cellindex_ty
if ( index == -1 ) break;
if ( kci == 1 )
{
index = cellindex[0]->data[index];
*sum += array[index];
*sumq += array[index]*array[index];
*n += 1;
......@@ -256,10 +256,13 @@ void compute_sum(int i, int *n, double *sum, double *sumq, int kci, cellindex_ty
static
void samplegrid(double missval, int nci, cellindex_type **cellindex, double *array1, double *array2, double *array3)
{
static bool lstat = true;
int kci = nci-1;
int ncells2 = cellindex[kci]->ncells;
int nx = 0;
double x = 0;
#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(missval, ncells2, kci, cellindex, array1, array2, array3)
//#pragma omp parallel for default(none) shared(missval, ncells2, kci, cellindex, array1, array2, array3)
#endif
for ( int i = 0; i < ncells2; ++i )
{
......@@ -270,7 +273,9 @@ void samplegrid(double missval, int nci, cellindex_type **cellindex, double *arr
double var1 = (n*n > n) ? (sumq*n - sum*sum) / (n*n - n) : missval;
if ( var1 < 0 && var1 > -1.e-5 ) var1 = 0;
array3[i] = var_to_std(var1, missval); // std1
if ( lstat && n ) { nx++; x+=n; }
}
if ( cdoVerbose && lstat ) { lstat = false; cdoPrint("Mean number of childs %g", nx?x/nx:0); }
}
......@@ -306,6 +311,14 @@ void *Samplegridicon(void *argument)
int streamID1 = streamOpenRead(cdoStreamName(0));
int vlistID1 = streamInqVlist(streamID1);
int gridsize = vlistGridsizeMax(vlistID1);
if ( cdoVerbose ) cdoPrint("Source gridsize = %d", gridsize);
if ( gridsize != cellindex[0]->ncells )
cdoAbort("Gridsize (%d) of input stream and first grid (%d) differ!", gridsize, cellindex[0]->ncells);
if ( vlistNumber(vlistID1) != CDI_REAL ) gridsize *= 2;
double *array1 = (double *) Malloc(gridsize*sizeof(double));
int vlistID2 = vlistDuplicate(vlistID1);
int vlistID3 = vlistDuplicate(vlistID1);
......@@ -333,10 +346,6 @@ void *Samplegridicon(void *argument)
int streamID3 = streamOpenWrite(cdoStreamName(2), cdoFiletype());
streamDefVlist(streamID3, vlistID3);
int gridsize = vlistGridsizeMax(vlistID1);
if ( vlistNumber(vlistID1) != CDI_REAL ) gridsize *= 2;
double *array1 = (double *) Malloc(gridsize*sizeof(double));
int gridsize2 = gridInqSize(gridID2);
if ( cdoVerbose ) cdoPrint("Target gridsize = %d", gridsize2);
if ( vlistNumber(vlistID2) != CDI_REAL ) gridsize2 *= 2;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment