From a47d36b4ec82c8bf59fbedb4ce2571f454478394 Mon Sep 17 00:00:00 2001 From: Tom Quinn Date: Thu, 13 Mar 2025 14:08:16 -0700 Subject: [PATCH] Added 2D Ewald for patch simulations. --- Makefile | 7 +++++-- ewald.h | 3 +++ kd.c | 64 +++++++++++++++++++++++++++++++++++++++++++++++++++----- kd.h | 7 +++++-- main.c | 49 ++++++++++++++++++++++++++++++++++--------- 5 files changed, 111 insertions(+), 19 deletions(-) diff --git a/Makefile b/Makefile index d4a55c5..2e40e34 100644 --- a/Makefile +++ b/Makefile @@ -9,8 +9,8 @@ default: direct clean: rm -f *.o -direct: main.o kd.o grav.o ewald.o - $(CC) $(CFLAGS) -o direct main.o kd.o grav.o ewald.o $(LIBS) +direct: main.o kd.o grav.o ewald.o ewald2.o + $(CC) $(CFLAGS) -o direct main.o kd.o grav.o ewald.o ewald2.o $(LIBS) main.o: main.c kd.h $(CC) $(CFLAGS) -c main.c @@ -24,3 +24,6 @@ grav.o: grav.c grav.h kd.h spline.h ewald.o: ewald.c ewald.h kd.h spline.h $(CC) $(CFLAGS) -c ewald.c +ewald2.o: ewald2.c ewald.h kd.h spline.h + $(CC) $(CFLAGS) -c ewald2.c + diff --git a/ewald.h b/ewald.h index adaf9c9..638bf73 100644 --- a/ewald.h +++ b/ewald.h @@ -7,4 +7,7 @@ void diaEwald(PARTICLE *,int,int,double); void blkEwald(PARTICLE *,int,PARTICLE *,int,int,double); void umkEwald(PARTICLE *,int,PARTICLE *,int,int,double); +void diaEwald2D(PARTICLE *,int,int,float *); +void blkEwald2D(PARTICLE *,int,PARTICLE *,int,int,float *); +void umkEwald2D(PARTICLE *,int,PARTICLE *,int,int,float *); #endif diff --git a/kd.c b/kd.c index 01c8953..186c369 100644 --- a/kd.c +++ b/kd.c @@ -255,7 +255,7 @@ void kdOrder(KD kd) /* ** A simple driver. No support for check-pointing here! */ -void kdGravSimple(KD kd,int iSoftType,int bPeriodic) +void kdGravSimple(KD kd,int iSoftType,int bPeriodic, int bPeriodic2) { PARTICLE *p; int i,n; @@ -274,6 +274,9 @@ void kdGravSimple(KD kd,int iSoftType,int bPeriodic) if (bPeriodic) { diaEwald(p,n,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + diaEwald2D(p,n,iSoftType,kd->fPeriod); + } else { diaGrav(p,n,iSoftType); } @@ -309,6 +312,7 @@ void WriteChkpt(KD kd,int lBlock,int kBlock) } h.iSoftType = kd->iSoftType; h.bPeriodic = kd->bPeriodic; + h.bPeriodic2 = kd->bPeriodic2; h.iBlockSize = kd->iBlockSize; /* ** Store time used up to now. @@ -361,6 +365,7 @@ void ReadChkpt(KD kd,int *plBlock,int *pkBlock) } kd->iSoftType = h.iSoftType; kd->bPeriodic = h.bPeriodic; + kd->bPeriodic2 = h.bPeriodic2; kd->iBlockSize = h.iBlockSize; /* ** Fix up CPU time. @@ -388,7 +393,7 @@ void ReadChkpt(KD kd,int *plBlock,int *pkBlock) ** Driver routine for direct. With support for check points! */ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, - int bVerbose) + int bPeriodic2, int bVerbose) { PARTICLE *p,*q; int k,l,n,nk,nb,nbk,bs,rs,rsk,iStride,bc,i; @@ -403,7 +408,7 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, nbk = nk/bs; rsk = nk%bs; iStride = 16384/bs; - if (bPeriodic) iStride = iStride/100; + if (bPeriodic || bPeriodic2) iStride = iStride/100; if (!iStride) iStride = 1; if (lBlock || kBlock) { l = lBlock; @@ -436,6 +441,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { diaEwald(&p[k*bs],bs,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + diaEwald2D(&p[k*bs],bs,iSoftType,kd->fPeriod); + } else { diaGrav(&p[k*bs],bs,iSoftType); } @@ -447,6 +455,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { diaEwald(&p[k*bs],rs,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + diaEwald2D(&p[k*bs],rs,iSoftType,kd->fPeriod); + } else { diaGrav(&p[k*bs],rs,iSoftType); } @@ -474,6 +485,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { blkEwald(&p[(l-1)*bs],bs,&p[k*bs],bs,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + blkEwald2D(&p[(l-1)*bs],bs,&p[k*bs],bs,iSoftType,kd->fPeriod); + } else { blkGrav(&p[(l-1)*bs],bs,&p[k*bs],bs,iSoftType); } @@ -485,6 +499,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { blkEwald(&p[l*bs],bs,&p[nb*bs],rs,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + blkEwald2D(&p[l*bs],bs,&p[nb*bs],rs,iSoftType,kd->fPeriod); + } else { blkGrav(&p[l*bs],bs,&p[nb*bs],rs,iSoftType); } @@ -519,6 +536,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { umkEwald(&p[l*bs],bs,&q[k*bs],bs,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + umkEwald2D(&p[l*bs],bs,&q[k*bs],bs,iSoftType,kd->fPeriod); + } else { umkGrav(&p[l*bs],bs,&q[k*bs],bs,iSoftType); } @@ -526,6 +546,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { umkEwald(&p[l*bs],bs,&q[k*bs],rsk,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + umkEwald2D(&p[l*bs],bs,&q[k*bs],rsk,iSoftType,kd->fPeriod); + } else { umkGrav(&p[l*bs],bs,&q[k*bs],rsk,iSoftType); } @@ -534,6 +557,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { umkEwald(&p[l*bs],rs,&q[k*bs],bs,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic2) { + umkEwald2D(&p[l*bs],rs,&q[k*bs],bs,iSoftType,kd->fPeriod); + } else { umkGrav(&p[l*bs],rs,&q[k*bs],bs,iSoftType); } @@ -541,6 +567,9 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, if (bPeriodic) { umkEwald(&p[l*bs],rs,&q[k*bs],rsk,iSoftType,kd->fPeriod[0]); } + else if (bPeriodic) { + umkEwald2D(&p[l*bs],rs,&q[k*bs],rsk,iSoftType,kd->fPeriod); + } else { umkGrav(&p[l*bs],rs,&q[k*bs],rsk,iSoftType); } @@ -555,12 +584,14 @@ void Grav(KD kd,int lBlock,int kBlock,int iSoftType,int bPeriodic, void kdGrav(KD kd,int iBlockSize,int iSoftType,int bPeriodic, + int bPeriodic2, int bVerbose) { kd->iSoftType = iSoftType; kd->bPeriodic = bPeriodic; + kd->bPeriodic2 = bPeriodic2; kd->iBlockSize = iBlockSize; - Grav(kd,0,0,iSoftType,bPeriodic,bVerbose); + Grav(kd,0,0,iSoftType,bPeriodic,bPeriodic2, bVerbose); } @@ -569,7 +600,7 @@ void kdRestart(KD kd,int bVerbose) int l,k; ReadChkpt(kd,&l,&k); - Grav(kd,l,k,kd->iSoftType,kd->bPeriodic,bVerbose); + Grav(kd,l,k,kd->iSoftType,kd->bPeriodic,kd->bPeriodic2,bVerbose); } @@ -623,6 +654,29 @@ void kdOutAccel(KD kd,char *pszFile) fclose(fp); } +void kdOutAccel1D(KD kd,char *pszFile, int iDir) +{ + FILE *fp; + int i,iCnt; + + fp = fopen(pszFile,"w"); + assert(fp != NULL); + fprintf(fp,"%d\n",kd->nParticles); + iCnt = 0; + for (i=0;inGas;++i) { + if (kd->bGas) fprintf(fp,"%.17g\n",kd->p[iCnt++].a[iDir]); + else fprintf(fp,"0\n"); + } + for (i=0;inDark;++i) { + if (kd->bDark) fprintf(fp,"%.17g\n",kd->p[iCnt++].a[iDir]); + else fprintf(fp,"0\n"); + } + for (i=0;inStar;++i) { + if (kd->bStar) fprintf(fp,"%.17g\n",kd->p[iCnt++].a[iDir]); + else fprintf(fp,"0\n"); + } + fclose(fp); + } void kdOutPot(KD kd,char *pszFile) { diff --git a/kd.h b/kd.h index 5cbbb02..a098009 100644 --- a/kd.h +++ b/kd.h @@ -39,6 +39,7 @@ typedef struct kdContext { char *pszChkptName; int iSoftType; int bPeriodic; + int bPeriodic2; int iBlockSize; int uSecond; int uMicro; @@ -61,6 +62,7 @@ struct chkptHeader { float fCenter[3]; int iSoftType; int bPeriodic; + int bPeriodic2; int iBlockSize; int uSecond; int uMicro; @@ -76,10 +78,11 @@ void kdSetSoft(KD,float); void kdInMark(KD,char *); void kdMarkOrder(KD); void kdOrder(KD); -void kdGravSimple(KD,int,int); -void kdGrav(KD,int,int,int,int); +void kdGravSimple(KD,int,int,int); +void kdGrav(KD,int,int,int,int,int); void kdRestart(KD,int); void kdOutAccel(KD,char *); +void kdOutAccel1D(KD kd,char *pszFile, int iDir); void kdOutPot(KD,char *); void kdFinish(KD); diff --git a/main.c b/main.c index 3843266..c91a0e3 100644 --- a/main.c +++ b/main.c @@ -12,6 +12,7 @@ void usage(void) fprintf(stderr,"USAGE:\n"); fprintf(stderr,"direct [-e ] [-uniform] [-plummer] [-spline]\n"); fprintf(stderr," [-p ] [-G ] [-b ]\n"); + fprintf(stderr," [-p2 ] [-splitout] \n"); fprintf(stderr," [-dgs] [-do ] [-o ] [-std]\n"); fprintf(stderr," [-i ] [-restart] [-v] [-t]\n"); fprintf(stderr,"Reads tipsy binary format from stdin\n"); @@ -19,7 +20,7 @@ void usage(void) exit(1); } -void main(int argc,char **argv) +int main(int argc,char **argv) { KD kd; int i,iBlockSize,bVerbose,bSoft,bTestGrav; @@ -27,6 +28,8 @@ void main(int argc,char **argv) int bStandard; /* Use tipsy standard format */ int iSoftType,bMark; int j,bPeriodic; + int bPeriodic2; /* 2 D periodicity */ + int bSplitOut; /* Split accelerations into 3 files */ int sec,usec; double G; float fSoft,fPeriod[3],fCenter[3]; @@ -47,6 +50,7 @@ void main(int argc,char **argv) G = 1.0; iBlockSize = 256; bPeriodic = 0; + bPeriodic2 = 0; for (j=0;j<3;++j) { fCenter[j] = 0.0; fPeriod[j] = FLT_MAX; @@ -104,6 +108,13 @@ void main(int argc,char **argv) fPeriod[2] = atof(argv[i]); ++i; } + else if (!strcmp(argv[i],"-p2")) { + ++i; + bPeriodic2 = 1; + fPeriod[0] = atof(argv[i]); + fPeriod[1] = atof(argv[i]); + ++i; + } else if (!strcmp(argv[i],"-do")) { ++i; if (i >= argc) usage(); @@ -131,6 +142,10 @@ void main(int argc,char **argv) bStandard = 1; ++i; } + else if (!strcmp(argv[i],"-splitout")) { + bSplitOut = 1; + ++i; + } else if (*argv[i] == '-') { p = argv[i]; ++p; @@ -160,6 +175,7 @@ void main(int argc,char **argv) } else usage(); } + assert(!bPeriodic || !bPeriodic2); if (bVerbose) { printf("DIRECT v1.2: Joachim Stadel, Jan 1995\n"); fflush(stdout); @@ -196,11 +212,11 @@ void main(int argc,char **argv) } kdTime(kd,&sec,&usec); if (bTestGrav) { - kdGravSimple(kd,iSoftType,bPeriodic); - } + kdGravSimple(kd,iSoftType,bPeriodic,bPeriodic2); + } else { - kdGrav(kd,iBlockSize,iSoftType,bPeriodic,bVerbose); - } + kdGrav(kd,iBlockSize,iSoftType,bPeriodic,bPeriodic2,bVerbose); + } } if (bVerbose) { kdTime(kd,&sec,&usec); @@ -208,9 +224,22 @@ void main(int argc,char **argv) fflush(stdout); } if (bMark) kdOrder(kd); - strcpy(achFile,ach); - strcat(achFile,".acc"); - kdOutAccel(kd,achFile); + if(bSplitOut) { + strcpy(achFile,ach); + strcat(achFile,".accx"); + kdOutAccel1D(kd,achFile, 0); + strcpy(achFile,ach); + strcat(achFile,".accy"); + kdOutAccel1D(kd,achFile, 1); + strcpy(achFile,ach); + strcat(achFile,".accz"); + kdOutAccel1D(kd,achFile, 2); + } + else { + strcpy(achFile,ach); + strcat(achFile,".acc"); + kdOutAccel(kd,achFile); + } strcpy(achFile,ach); strcat(achFile,".pot"); kdOutPot(kd,achFile); @@ -227,7 +256,7 @@ void main(int argc,char **argv) } } kdFinish(kd); - } - + return 0; + }