mirror of
https://bitbucket.org/cosmicvoids/vide_public.git
synced 2025-07-04 15:21:11 +00:00
Harden zobov for problems link to single precisions (not fool proof yet but make the problem more unlikely)
This commit is contained in:
Guilhem Lavaux
parent
f1e80ac6ed
commit
aaef7634cd
3 changed files with 72 additions and 82 deletions
|
|
@ -124,19 +124,26 @@ int main(int argc,char **argv) {
|
|||
fread(&j,1,sizeof(pid_t),adj);
|
||||
if (j < np) {
|
||||
/* Set both halves of the pair */
|
||||
assert(i < j);
|
||||
assert(i < j);
|
||||
if (p[i].ncnt == p[i].nadj)
|
||||
{
|
||||
p[i].adj = (pid_t *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(pid_t));
|
||||
p[i].nadj++;
|
||||
int q;
|
||||
printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", i, j);
|
||||
for (q=0;q<p[i].nadj;q++)
|
||||
printf(" %d\n", p[i].adj[q]);
|
||||
abort();
|
||||
}
|
||||
if (p[j].ncnt == p[j].nadj)
|
||||
if (p[j].ncnt == p[j].nadj)
|
||||
{
|
||||
p[j].adj = (pid_t *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(pid_t));
|
||||
p[j].nadj++;
|
||||
}
|
||||
p[i].adj[p[i].ncnt] = j;
|
||||
p[j].adj[p[j].ncnt] = i;
|
||||
int q;
|
||||
printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", j, i);
|
||||
for (q=0;q<p[j].nadj;q++)
|
||||
printf(" %d\n", p[j].adj[q]);
|
||||
abort();
|
||||
}
|
||||
|
||||
p[i].adj[p[i].ncnt] = j;
|
||||
p[j].adj[p[j].ncnt] = i;
|
||||
p[i].ncnt++; p[j].ncnt++;
|
||||
} else {
|
||||
printf("%d: adj = %d\n",i,j);
|
||||
|
|
|
|||
|
|
@ -21,16 +21,16 @@ int main(int argc, char *argv[]) {
|
|||
char *posfile, outfile[200], *suffix, *outDir;
|
||||
PARTADJ *adjs;
|
||||
float *vols;
|
||||
float predict, xmin,xmax,ymin,ymax,zmin,zmax;
|
||||
realT predict, xmin,xmax,ymin,ymax,zmin,zmax;
|
||||
pid_t *orig;
|
||||
|
||||
int isitinbuf;
|
||||
char isitinmain, d;
|
||||
int numdiv;
|
||||
pid_t nvp, nvpall, nvpbuf;
|
||||
float width, width2, totwidth, totwidth2, bf, s, g;
|
||||
realT width, width2, totwidth, totwidth2, bf, s, g;
|
||||
float border, boxsize;
|
||||
float c[3];
|
||||
realT c[3];
|
||||
int b[3];
|
||||
double totalvol;
|
||||
|
||||
|
|
@ -118,7 +118,7 @@ int main(int argc, char *argv[]) {
|
|||
exit(0);
|
||||
}
|
||||
|
||||
DL c[d] = ((float)b[d]+0.5)*width;
|
||||
DL c[d] = ((float)b[d])*width;
|
||||
printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
|
||||
/* Assign temporary array*/
|
||||
nvpbuf = 0; /* Number of particles to tesselate, including
|
||||
|
|
@ -159,7 +159,7 @@ int main(int argc, char *argv[]) {
|
|||
for (i=0; i<np; i++) {
|
||||
isitinmain = 1;
|
||||
DL {
|
||||
rtemp[d] = r[i][d] - c[d];
|
||||
rtemp[d] = (realT)r[i][d] - (realT)c[d];
|
||||
if (rtemp[d] > 0.5) rtemp[d] --;
|
||||
if (rtemp[d] < -0.5) rtemp[d] ++;
|
||||
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
|
||||
|
|
@ -186,7 +186,7 @@ int main(int argc, char *argv[]) {
|
|||
isitinmain = 1;
|
||||
|
||||
DL {
|
||||
rtemp[d] = r[i][d] - c[d];
|
||||
rtemp[d] = (realT)r[i][d] - (realT)c[d];
|
||||
if (rtemp[d] > 0.5) rtemp[d] --;
|
||||
if (rtemp[d] < -0.5) rtemp[d] ++;
|
||||
isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2);
|
||||
|
|
@ -223,40 +223,40 @@ int main(int argc, char *argv[]) {
|
|||
for (i=0; i<NGUARD+1; i++) {
|
||||
for (j=0; j<NGUARD+1; j++) {
|
||||
/* Bottom */
|
||||
parts[3*nvpall] = -width2 + (float)i * s;
|
||||
parts[3*nvpall+1] = -width2 + (float)j * s;
|
||||
parts[3*nvpall] = -width2 + (realT)i * s;
|
||||
parts[3*nvpall+1] = -width2 + (realT)j * s;
|
||||
parts[3*nvpall+2] = -width2 - g;
|
||||
nvpall++;
|
||||
/* Top */
|
||||
parts[3*nvpall] = -width2 + (float)i * s;
|
||||
parts[3*nvpall+1] = -width2 + (float)j * s;
|
||||
parts[3*nvpall] = -width2 + (realT)i * s;
|
||||
parts[3*nvpall+1] = -width2 + (realT)j * s;
|
||||
parts[3*nvpall+2] = width2 + g;
|
||||
nvpall++;
|
||||
}
|
||||
}
|
||||
for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
|
||||
for (j=0; j<NGUARD+1; j++) {
|
||||
parts[3*nvpall] = -width2 + (float)i * s;
|
||||
parts[3*nvpall] = -width2 + (realT)i * s;
|
||||
parts[3*nvpall+1] = -width2 - g;
|
||||
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||
parts[3*nvpall+2] = -width2 + (realT)j * s;
|
||||
nvpall++;
|
||||
|
||||
parts[3*nvpall] = -width2 + (float)i * s;
|
||||
parts[3*nvpall] = -width2 + (realT)i * s;
|
||||
parts[3*nvpall+1] = width2 + g;
|
||||
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||
parts[3*nvpall+2] = -width2 + (realT)j * s;
|
||||
nvpall++;
|
||||
}
|
||||
}
|
||||
for (i=0; i<NGUARD+1; i++) {
|
||||
for (j=0; j<NGUARD+1; j++) {
|
||||
parts[3*nvpall] = -width2 - g;
|
||||
parts[3*nvpall+1] = -width2 + (float)i * s;
|
||||
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||
parts[3*nvpall+1] = -width2 + (realT)i * s;
|
||||
parts[3*nvpall+2] = -width2 + (realT)j * s;
|
||||
nvpall++;
|
||||
|
||||
parts[3*nvpall] = width2 + g;
|
||||
parts[3*nvpall+1] = -width2 + (float)i * s;
|
||||
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||
parts[3*nvpall+1] = -width2 + (realT)i * s;
|
||||
parts[3*nvpall+2] = -width2 + (realT)j * s;
|
||||
nvpall++;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,3 +1,4 @@
|
|||
#include <assert.h>
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
|
@ -96,7 +97,11 @@ int main(int argc, char *argv[]) {
|
|||
exit(0);
|
||||
}
|
||||
for (p=0;p<np;p++)
|
||||
vols[p] = -1.;
|
||||
{
|
||||
vols[p] = -1.;
|
||||
adjs[p].nadj = 0;
|
||||
adjs[p].adj = 0;
|
||||
}
|
||||
|
||||
nvpsum = 0;
|
||||
for (i = 0; i < numdiv; i++) {
|
||||
|
|
@ -118,13 +123,9 @@ int main(int argc, char *argv[]) {
|
|||
for (p=0;p<nvp;p++) {
|
||||
fread(&volstemp,1,sizeof(float),part);
|
||||
if (vols[orig[p]] > -1.)
|
||||
// PMS
|
||||
if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3 && orig[p] < mockIndex) {
|
||||
//if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) {
|
||||
// END PMS
|
||||
printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n",
|
||||
orig[p],vols[orig[p]],volstemp);
|
||||
// TEST
|
||||
//exit(0);
|
||||
}
|
||||
vols[orig[p]] = volstemp;
|
||||
|
|
@ -132,10 +133,13 @@ int main(int argc, char *argv[]) {
|
|||
|
||||
for (p=0;p<nvp;p++) {
|
||||
fread(&na,1,sizeof(int),part);
|
||||
pid_t pid = orig[p];
|
||||
if (na > 0) {
|
||||
assert(adjs[pid].nadj == 0);// || adjs[pid].nadj == na);
|
||||
adjs[orig[p]].nadj = na;
|
||||
adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t));
|
||||
if (adjs[orig[p]].adj == NULL) {
|
||||
if (adjs[pid].adj == 0)
|
||||
adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t));
|
||||
if (adjs[orig[p]].adj == 0) {
|
||||
printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
|
||||
exit(0);
|
||||
}
|
||||
|
|
@ -160,51 +164,37 @@ int main(int argc, char *argv[]) {
|
|||
adjs[i].nadj = 0;
|
||||
}
|
||||
|
||||
// unlink particles adjacent to mock galaxies
|
||||
for (i = 0; i < mockIndex; i++) {
|
||||
for (j = 0; j < adjs[i].nadj; j++) {
|
||||
if (adjs[i].adj[j] > mockIndex) {
|
||||
//printf("KILLING %d\n", i);
|
||||
vols[i] = 1.e-29;
|
||||
adjs[i].nadj = 0;
|
||||
numRemoved++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// update all other adjacencies
|
||||
for (i = 0; i < mockIndex; i++) {
|
||||
|
||||
int numAdjSaved = 0;
|
||||
for (j = 0; j < adjs[i].nadj; j++) {
|
||||
|
||||
//if ( vols[adjs[i].adj[j]] != -99) {
|
||||
if ( adjs[adjs[i].adj[j]].nadj != 0) {
|
||||
adjs[i].adj[numAdjSaved] = adjs[i].adj[j];
|
||||
numAdjSaved++;
|
||||
}
|
||||
|
||||
}
|
||||
adjs[i].nadj = numAdjSaved;
|
||||
|
||||
}
|
||||
|
||||
/*
|
||||
// unlink particles adjacent to mock galaxies
|
||||
for (i = 0; i < mockIndex; i++) {
|
||||
printf("ADJ: %d %d : ", i, adjs[i].nadj);
|
||||
for (j = 0; j < adjs[i].nadj; j++) {
|
||||
printf(" %d", adjs[i].adj[j]);
|
||||
if (adjs[i].adj[j] > mockIndex) {
|
||||
//printf("KILLING %d\n", i);
|
||||
vols[i] = 1.e-29;
|
||||
adjs[i].nadj = 0;
|
||||
numRemoved++;
|
||||
break;
|
||||
}
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
*/
|
||||
|
||||
|
||||
// update all other adjacencies
|
||||
for (i = 0; i < mockIndex; i++) {
|
||||
|
||||
int numAdjSaved = 0;
|
||||
for (j = 0; j < adjs[i].nadj; j++) {
|
||||
|
||||
printf("Removed %d mock galaxies and %d adjacent galaxies.\n", np-mockIndex,
|
||||
if ( adjs[adjs[i].adj[j]].nadj != 0) {
|
||||
adjs[i].adj[numAdjSaved] = adjs[i].adj[j];
|
||||
numAdjSaved++;
|
||||
}
|
||||
|
||||
}
|
||||
adjs[i].nadj = numAdjSaved;
|
||||
}
|
||||
|
||||
printf("Removed %d mock galaxies and %d adjacent galaxies.\n", np-mockIndex,
|
||||
numRemoved);
|
||||
printf("There are %d galaxies remaining.\n", mockIndex-numRemoved);
|
||||
printf("There are %d galaxies remaining.\n", mockIndex-numRemoved);
|
||||
|
||||
// END PMS
|
||||
|
||||
|
|
@ -238,10 +228,7 @@ printf("\n");
|
|||
printf("Unable to open %s\n",adjfile);
|
||||
exit(0);
|
||||
}
|
||||
// PMS
|
||||
fwrite(&mockIndex,1, sizeof(int),adj);
|
||||
//fwrite(&np,1, sizeof(int),adj);
|
||||
// END OMS
|
||||
/* Adjacencies: first the numbers of adjacencies,
|
||||
and the number we're actually going to write per particle */
|
||||
|
||||
|
|
@ -249,19 +236,15 @@ printf("\n");
|
|||
for(i=0;i<mockIndex;i++)
|
||||
cnt_adj[i] = adjs[i].nadj;
|
||||
|
||||
// PMS
|
||||
for (i=0;i<mockIndex;i++)
|
||||
// END PMS
|
||||
fwrite(&cnt_adj[i],1,sizeof(int),adj);
|
||||
|
||||
/* Now the lists of adjacencies (without double counting) */
|
||||
// PMS
|
||||
for (i=0;i<mockIndex;i++) {
|
||||
//for (i=0;i<np;i++) {
|
||||
//if (adjs[i].nadj > 0) {
|
||||
// END PMS
|
||||
nout = 0;
|
||||
for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++;
|
||||
for (j=0;j<adjs[i].nadj; j++)
|
||||
if (adjs[i].adj[j] > i)
|
||||
nout++;
|
||||
fwrite(&nout,1,sizeof(int),adj);
|
||||
for (j=0;j<adjs[i].nadj; j++) {
|
||||
pid_t id = adjs[i].adj[j];
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue