From 23d665f7bd9b0c594c02f198d8a7fa2f159ce1dd Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Thu, 22 May 2025 09:57:08 -0400 Subject: [PATCH] FIXED zobov bug when only using one subdomain; added zobov buffer as tunable parameter --- c_source/zobov/voz1b1.c | 76 ++++++++++----- c_source/zobov/vozinit.c | 119 +++++++++++++----------- c_source/zobov/voztie.c | 12 ++- python_source/backend/launchers.py | 10 +- python_source/vide_pipeline/__main__.py | 1 + python_source/voidUtil/catalogUtil.py | 4 +- 6 files changed, 139 insertions(+), 83 deletions(-) diff --git a/c_source/zobov/voz1b1.c b/c_source/zobov/voz1b1.c index 23a78ab..62e5a0a 100644 --- a/c_source/zobov/voz1b1.c +++ b/c_source/zobov/voz1b1.c @@ -17,7 +17,7 @@ int main(int argc, char *argv[]) { coordT rtemp[3], *parts; coordT deladjs[3*MAXVERVER], points[3*MAXVERVER]; pointT intpoints[3*MAXVERVER]; - FILE *pos, *out; + FILE *pos, *out, *testFile; char *posfile, outfile[200], *suffix, *outDir; PARTADJ *adjs; float *vols; @@ -33,8 +33,11 @@ int main(int argc, char *argv[]) { realT c[3]; int b[3]; double totalvol; + int isObs = 0; - if (argc != 10) { + printf("Routine: voz1b1\n"); + + if (argc != 11) { printf("Wrong number of arguments.\n"); printf("arg1: position file\n"); printf("arg2: border size\n"); @@ -43,6 +46,7 @@ int main(int argc, char *argv[]) { printf("arg5: number of divisions\n"); printf("arg6-8: b[0-2]\n\n"); printf("arg9: output directory\n"); + printf("arg10: observation flag\n"); exit(0); } posfile = argv[1]; @@ -79,8 +83,15 @@ int main(int argc, char *argv[]) { exit(0); } outDir = argv[9]; - + if (sscanf(argv[10],"%d",&isObs) != 1) { + printf("That's no observation mode; try again.\n"); + exit(0); + } + + printf("Working with subdomain b=[%d,%d,%d]\n", b[0], b[1], b[2]); + /* Boxsize should be the range in r, yielding a range 0-1 */ + printf("Reading particle file...\n"); np = posread(posfile,&r,1./boxsize); printf("%d particles\n",np);fflush(stdout); @@ -90,7 +101,8 @@ int main(int argc, char *argv[]) { if (r[i][1]ymax) ymax = r[i][1]; if (r[i][2]zmax) zmax = r[i][2]; } - printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout); + printf("Full box particle stats: np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np, + xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout); width = 1./(float)numdiv; width2 = 0.5*width; @@ -108,7 +120,7 @@ int main(int argc, char *argv[]) { exit(0); } g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf))); - printf("s = %f, bf = %f, g = %f.\n",s,bf,g); + printf("Guard particle information: s = %f, bf = %f, g = %f.\n",s,bf,g); fflush(stdout); @@ -118,11 +130,10 @@ int main(int argc, char *argv[]) { exit(1); } - DL c[d] = ((float)b[d])*width; - printf("c: %f,%f,%f\n",c[0],c[1],c[2]); + DL c[d] = ((float)b[d]+0.5)*width; + printf("Counting particles in subdomain...\n"); /* Assign temporary array*/ - nvpbuf = 0; /* Number of particles to tesselate, including - buffer */ + nvpbuf = 0; /* Number of particles to tesselate, includng ibuffer */ nvp = 0; /* Without the buffer */ for (i=0; i zmax) zmax = rtemp[2]; } } - printf("nvp = %d\n",nvp); - printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax); + printf("Subdomain: b=(%d,%d,%d), c=(%f,%f,%f)\n", + b[0],b[1],b[2],c[0],c[1],c[2]); + printf(" nvp=%d, nvpbuf=%d\n", nvp,nvpbuf); + printf(" Particle ranges: x: %f,%f; y: %f,%f; z:%f,%f\n", + xmin,xmax,ymin,ymax,zmin,zmax); + + printf("Adding particles to buffer regions...\n"); nvpbuf = nvp; for (i=0; i zmax) zmax = rtemp[2]; } } - printf("nvpbuf = %d\n",nvpbuf); - printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax); nvpall = nvpbuf; predict = pow(width+2.*bf,3)*(float)np; - printf("There should be ~ %f points; there are %d\n",predict,nvpbuf); + printf("Total particles including buffer = %d (predicted ~%f)\n", + nvpbuf, predict); + printf(" Particle ranges including buffer: x: %f,%f; y: %f,%f; z:%f,%f\n", + xmin,xmax,ymin,ymax,zmin,zmax); for (i=0;i nvp) itype = 1; // it's a buffer particle + if (i > nvpbuf) itype = 2; // it's a guard particle + fprintf(testFile, "%f %f %f %d %f %f %f\n", c[0], c[1], c[2], + itype, + parts[3*i], + parts[3*i+1], + parts[3*i+2]); + } + fclose(testFile); +*/ + xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; for (i=nvpbuf;i 0.5) rtemp[d] --; - if (rtemp[d] < -0.5) rtemp[d] ++; - isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2); - isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); - } - if (isitinbuf) { - nvpbuf++; - } - if (isitinmain) { - nvp++; - if (rtemp[0] < xmin) xmin = rtemp[0]; - if (rtemp[0] > xmax) xmax = rtemp[0]; - if (rtemp[1] < ymin) ymin = rtemp[1]; - if (rtemp[1] > ymax) ymax = rtemp[1]; - if (rtemp[2] < zmin) zmin = rtemp[2]; - if (rtemp[2] > zmax) zmax = rtemp[2]; - } - } - if (nvp > nvpmax) nvpmax = nvp; - if (nvpbuf > nvpbufmax) nvpbufmax = nvpbuf; - if (nvp < nvpmin) nvpmin = nvp; - if (nvpbuf < nvpbufmin) nvpbufmin = nvpbuf; - - printf("b=(%d,%d,%d), c=(%f,%f,%f), nvp=%d, nvpbuf=%d\n", - b[0],b[1],b[2],c[0],c[1],c[2],nvp,nvpbuf); + nvp = 0; /* Number of particles excluding buffer */ + nvpbuf = 0; /* Number of particles to tesselate, including buffer */ + xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; + + for (i=0; i 0.5) rtemp[d] --; + if (rtemp[d] < -0.5) rtemp[d] ++; + //} + isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2); + isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); + } + if (isitinbuf) { nvpbuf++; } + if (isitinmain) { + nvp++; + if (rtemp[0] < xmin) xmin = rtemp[0]; + if (rtemp[0] > xmax) xmax = rtemp[0]; + if (rtemp[1] < ymin) ymin = rtemp[1]; + if (rtemp[1] > ymax) ymax = rtemp[1]; + if (rtemp[2] < zmin) zmin = rtemp[2]; + if (rtemp[2] > zmax) zmax = rtemp[2]; + } } - } + + if (nvp > nvpmax) nvpmax = nvp; + if (nvpbuf > nvpbufmax) nvpbufmax = nvpbuf; + if (nvp < nvpmin) nvpmin = nvp; + if (nvpbuf < nvpbufmin) nvpbufmin = nvpbuf; + + printf("Subdomain: b=(%d,%d,%d), c=(%f,%f,%f)\n", + b[0],b[1],b[2],c[0],c[1],c[2]); + printf(" nvp=%d, nvpbuf=%d\n", nvp,nvpbuf); + printf(" Particle ranges: x: %f,%f; y: %f,%f; z:%f,%f\n", + xmin,xmax,ymin,ymax,zmin,zmax); + } + } } printf("Nvp range: %d,%d\n",nvpmin,nvpmax); printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax); numGuards = 6*(NGUARD+1)*(NGUARD+1); + printf("Num guards = %d\n", numGuards); printf("Total max particles: %d\n" , nvpbufmax+numGuards); - if (nvpbufmax+numGuards >= QHULL_MAX_PARTICLES) - { + if (nvpbufmax+numGuards >= QHULL_MAX_PARTICLES) { printf("Too many particles to tesselate per division (Qhull will overflow). Increase divisions or reduce buffer size.\n"); fflush(stdout); exit(1); - } + } /* Output script file */ sprintf(scrfile,"scr%s",suffix); @@ -164,23 +179,19 @@ int main(int argc, char *argv[]) { fprintf(scr,"#!/bin/bash -f\n"); p = 0; for (b[0]=0;b[0]