diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc index 86b8a7680efca451f6c3fb9d26a1620f07124684..897b1319c2828eea09bdf73956b0263aa6cd0ab4 100644 --- a/tracy/tools/soltracy.cc +++ b/tracy/tools/soltracy.cc @@ -9,7 +9,17 @@ int no_tps = ORDER; // arbitrary TPSA order is defined locally extern bool freq_map; -#include "tracy_lib.h" +#if HAVE_CONFIG_H +#include <config.h> +#endif + +#if MPI_EXEC + //library of parallel computation,must be before "stdio.h" + #include <mpi.h> +#endif + + #include "tracy_lib.h" + //*************************************************************************************** // @@ -32,6 +42,7 @@ int main(int argc, char *argv[]) { /* parameters to read the user input script .prm */ long i=0L; //initialize the for loop to read command string + long j=0L; //initialize the for loop to read the files from parallel computing char CommandStr[max_str]; double nux = 0.0, nuy = 0.0, ksix = 0.0, ksiy = 0.0; bool chroma=true; @@ -46,7 +57,10 @@ int main(int argc, char *argv[]) { const int NCOMMAND = 500; UserCommand UserCommandFlag[NCOMMAND]; - +#if MPI_EXEC + //Initialize parallel computation + MPI_Init(&argc, &argv); +#endif /************************************************************************ read in files and flags *************************************************************************/ @@ -126,6 +140,13 @@ int main(int argc, char *argv[]) { cout << "globval.quad_fringe is " << globval.quad_fringe << "\n"; } + //deactive quadrupole fringe field + else if(strcmp(CommandStr,"QuadFringeOffFlag") == 0) { + globval.quad_fringe = false; + cout << "\n"; + cout << "globval.quad_fringe is " << globval.quad_fringe << "\n"; + } + //set RF voltage else if(strcmp(CommandStr,"RFvoltageFlag") == 0) { printf("\nSetting RF voltage:\n"); @@ -141,133 +162,60 @@ int main(int argc, char *argv[]) { ReadCh(UserCommandFlag[i].chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/ PrintCh(); // print chamber into chamber.out } - - // read the misalignment errors to the elements, then do COD correction - // using SVD method. + + // read the misalignment errors to the elements // Based on the function error_and_correction() in nsls-ii_lib_templ.h - - // else if(strcmp(ReadaefileFlag, "") != 0){ else if (strcmp(CommandStr,"ReadaefileFlag") == 0) { - // ***** Apply corrections and output flatfile for n_stat sets of random #'s - bool cod = true; - int k, icod=0; - FILE *hOrbitFile, *vOrbitFile ; - int hcorrIdx[nCOR], vcorrIdx[nCOR]; //list of corr for orbit correction - - //initialize the corrector list - for ( k = 0; k < nCOR; k++){ - hcorrIdx[k] = -1; - vcorrIdx[k] = -1; - } - - //Get response matrix between bpm and correctors, and then print the SVD setting to the files - if (n_orbit!=0 || n_scale!=0) { - - // select correctors to be used - readCorrectorList(hcorr_file, vcorr_file, hcorrIdx, vcorrIdx); - - fprintf(stdout, "\n\nSVD correction setting:\n"); - fprintf(stdout, "H-plane %d singular values:\n", nwh); - fprintf(stdout, "V-plane %d singular values:\n\n",nwv); - - // compute beam response matrix - printf("\n"); - printf("Computing beam response matrix\n"); - //get the response matrix between bpm and correctors - gcmats(globval.bpm, globval.hcorr, 1, hcorrIdx); - gcmats(globval.bpm, globval.vcorr, 2, vcorrIdx); - /* gcmat(globval.bpm, globval.hcorr, 1); - gcmat(globval.bpm, globval.vcorr, 2);*/ - - // print response matrices to files 'svdh.out' and file 'svdv.out' - prt_gcmat(globval.bpm, globval.hcorr, 1); - prt_gcmat(globval.bpm, globval.vcorr, 2); - - } - - //print the statistics of orbit in file 'OrbScanFile.out' - OrbScanFile = file_write(OrbScanFileName); - - //write files with orbits at all element locations - hOrbitFile = file_write(hOrbitFileName); - vOrbitFile = file_write(vOrbitFileName); - - fprintf(hOrbitFile, "# First line: s-location (m) \n"); - fprintf(hOrbitFile, "# After orbit correction: Horizontal closed orbit at all element locations (with %3d BPMs) at different loop\n", GetnKid(globval.bpm)); - fprintf(vOrbitFile, "# First line s-location (m) \n"); - fprintf(vOrbitFile, "# After orbit correction: Vertical closed orbit at all element locations (with %3d BPMs) at different loop\n", GetnKid(globval.bpm)); - - for (k = 0; k < globval.Cell_nLoc; k++){ - fprintf(hOrbitFile, "% 9.3e ", Cell[k].S); - fprintf(vOrbitFile, "% 9.3e ", Cell[k].S); - } // end for - - fprintf(hOrbitFile, "\n"); - fprintf(vOrbitFile, "\n"); + // load misalignments + cout << "Read misalignment errors from file: " << UserCommandFlag[i].ae_file << endl; + LoadAlignTol(UserCommandFlag[i].ae_file, true, 1.0, true, 0); + Ring_GetTwiss(true, 0.0); + } + //do COD correction using SVD method. + else if(strcmp(CommandStr,"CODCorrectFlag") == 0) { + cout << "Begin Closed Orbit correction: " << endl; - for (k = 1; k <= n_stat; k++) {// loop for n_stat sets of errors - - if (n_orbit!=0 || n_scale!=0) { - - cod = - CorrectCOD_Ns(hOrbitFile, vOrbitFile, UserCommandFlag[i].ae_file, - n_orbit,n_scale,k,nwh,nwv, hcorrIdx, vcorrIdx); - /* cod = CorrectCOD_N(ae_file, n_orbit, n_scale, k);*/ - printf("\n"); - - if (cod){ - /* printf("done with orbit correction, now do coupling", - " correction plus vert. disp\n");*/ - if(n_orbit == 0) - printf("iter # %3d n_obit = 0, no orbit correction \n",k); - else - printf("iter # %3d Orbit correction succeeded\n", k); - } - else - if(!cod){ - icod = icod + 1; - fprintf(stdout, "!!! iter # %3d error_and_correction\n",k); - } - //chk_cod(cod, "iter # %3d error_and_correction"); - } - // Ring_GetTwiss(chroma, dp); - Ring_GetTwiss(true, 0.0); - - // for debugging - //print flat lattice - //sprintf(mfile_name, "flat_file.%03d.dat",k); - //prtmfile(mfile_name); - - } + if(n_orbit < 1){ + cout << "interation number n_obit should NOT small than 1!!!" << endl; + exit_(1); + } - fprintf(stdout, "Number of unstable orbits %d/%d", icod, n_stat); - prt_cod("corr_after.out", globval.bpm, true); - - // close file giving orbit at BPM location - fclose(hOrbitFile); - fclose(vOrbitFile); - fclose(OrbScanFile); + CODCorrect(hcorr_file,vcorr_file, n_orbit,nwh,nwv); + Ring_GetTwiss(true, 0.0); + prt_cod("summary_miserr_codcorr.out", globval.bpm, true); } - // set the field error into the lattice // The corresponding field error is replaced by the new value. // This feature is generic, works for all lattices - else if (strcmp(CommandStr,"ReadfefileFlag") == 0) { - fprintf(stdout,"\n Read field error from fe_file: %s \n", UserCommandFlag[i].fe_file); + else if (strcmp(CommandStr,"ReadfefileFlag") == 0) { + fprintf(stdout,"\n Read multipole field errors from file: %s \n", UserCommandFlag[i].fe_file); LoadFieldErrs(UserCommandFlag[i].fe_file, true, 1.0, true, 1); Ring_GetTwiss(true, 0.0); prtmfile("flat_file_fefile.dat"); } + // read multipole errors from a file; specific for soleil lattice else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) { - fprintf(stdout,"\n Read Multipoles file for lattice with thick sextupoles, specific for SOLEIL lattice: \n"); + fprintf(stdout,"\n Read Multipoles field error for SOLEIL lattice with thick sextupoles, from file:\n"); + fprintf(stdout," %s\n",multipole_file); ReadFieldErr(multipole_file); //first print the full lattice with error as a flat file prtmfile("flat_file_errmultipole.dat"); // writes flat file with all element errors /* very important file for debug*/ Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */ printglob(); } + // read the sources of coupling from a file; for soleil lattice + else if(strcmp(CommandStr,"ReadVirtualSkewquadFlag") == 0) { + fprintf(stdout,"\n Read virtual skew quadrupole setting of SOLEIL lattice from file:\n"); + fprintf(stdout," %s \n",virtualskewquad_file); + + ReadVirtualSkewQuad(virtualskewquad_file); + //first print the full lattice with sources of coupling as a flat file + prtmfile("flat_file_skewquad.dat"); /* very important file for debug*/ + //get the coupling + Coupling_Edwards_Teng(); + } //print the twiss paramaters in a file defined by the name else if(strcmp(CommandStr,"PrintTwissFlag") == 0) { @@ -431,8 +379,8 @@ int main(int argc, char *argv[]) { printlatt("linlat_coupling.out"); /* dump linear lattice functions into "linlat.dat" */ // GetEmittance(ElemIndex("cav"), true); //only for test Coupling_Edwards_Teng(); - Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */ - printglob(); /* print parameter list */ + // Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */ + // printglob(); /* print parameter list */ } // add coupling by random rotating of the full quadrupole magnets @@ -492,27 +440,162 @@ int main(int argc, char *argv[]) { // Computes FMA else if(strcmp(CommandStr,"FmapFlag") == 0) { printf("\n begin Fmap calculation for on momentum particles: \n"); - fmap(UserCommandFlag[i]._FmapFlag_fmap_file, - UserCommandFlag[i]._FmapFlag_nxpoint, - UserCommandFlag[i]._FmapFlag_nypoint, - UserCommandFlag[i]._FmapFlag_nturn, - UserCommandFlag[i]._FmapFlag_xmax, - UserCommandFlag[i]._FmapFlag_ymax, - UserCommandFlag[i]._FmapFlag_delta, - UserCommandFlag[i]._FmapFlag_diffusion); + + #if MPI_EXEC + + //initialization for parallel computing + int myid=0, numprocs=0; + int namelen=0; + char processor_name[MPI_MAX_PROCESSOR_NAME]=""; + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Get_processor_name(processor_name, &namelen); + printf("\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); + + printf("\n begin Fmap calculation for on momentum particles: \n"); + //Each core or process, characterized by myid, will track different region of fmap + fmap_p(UserCommandFlag[i]._FmapFlag_fmap_file, + UserCommandFlag[i]._FmapFlag_nxpoint, + UserCommandFlag[i]._FmapFlag_nypoint, + UserCommandFlag[i]._FmapFlag_nturn, + UserCommandFlag[i]._FmapFlag_xmax, + UserCommandFlag[i]._FmapFlag_ymax, + UserCommandFlag[i]._FmapFlag_delta, + UserCommandFlag[i]._FmapFlag_diffusion, + numprocs,myid); + + //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed. + MPI_Barrier(MPI_COMM_WORLD); + + //Collecting data from all files generated by cores, then write to fmap_p.out + if(myid==0) + { + // system("cp 0fmap.out fmap_p.out"); + + FILE *outf; + if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL) + { + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file); + exit_(1); + } + + FILE *fp; + char buffer[81]; + char FmapFile[50]; + for(int j=0; j<numprocs; j++) + { + FmapFile[0]='\0'; + sprintf(FmapFile,"%d",j); + strcat(FmapFile,UserCommandFlag[i]._FmapFlag_fmap_file); + + if((fp = fopen(FmapFile, "r")) == NULL) + { + fprintf(stdout, "%s: error while opening file.\n", FmapFile); + exit_(1); + } + + while(fgets(buffer, 80, fp) != NULL) + { + fputs(buffer, outf); + buffer[0]='\0'; + } + fclose(fp); + } + fclose(outf); + } + MPI_Barrier(MPI_COMM_WORLD); + printf("Process %d finished for on momentum fmap.\n", myid); + + #else + fmap(UserCommandFlag[i]._FmapFlag_fmap_file, + UserCommandFlag[i]._FmapFlag_nxpoint, + UserCommandFlag[i]._FmapFlag_nypoint, + UserCommandFlag[i]._FmapFlag_nturn, + UserCommandFlag[i]._FmapFlag_xmax, + UserCommandFlag[i]._FmapFlag_ymax, + UserCommandFlag[i]._FmapFlag_delta, + UserCommandFlag[i]._FmapFlag_diffusion); +#endif } // Compute FMA dp else if(strcmp(CommandStr,"FmapdpFlag") == 0) { printf("\n begin Fmap calculation for off momentum particles: \n"); - fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, - UserCommandFlag[i]._FmapdpFlag_nxpoint, - UserCommandFlag[i]._FmapdpFlag_nepoint, - UserCommandFlag[i]._FmapdpFlag_nturn, - UserCommandFlag[i]._FmapdpFlag_xmax, - UserCommandFlag[i]._FmapdpFlag_emax, - UserCommandFlag[i]._FmapdpFlag_z, - UserCommandFlag[i]._FmapdpFlag_diffusion); + + #if MPI_EXEC + + //initialization for parallel computing + int myid=0, numprocs=0; + int namelen=0; + char processor_name[MPI_MAX_PROCESSOR_NAME]=""; + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Get_processor_name(processor_name, &namelen); + printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); + + printf("\n begin Fmap calculation for off momentum particles: \n"); + fmapdp_p(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, + UserCommandFlag[i]._FmapdpFlag_nxpoint, + UserCommandFlag[i]._FmapdpFlag_nepoint, + UserCommandFlag[i]._FmapdpFlag_nturn, + UserCommandFlag[i]._FmapdpFlag_xmax, + UserCommandFlag[i]._FmapdpFlag_emax, + UserCommandFlag[i]._FmapdpFlag_z, + UserCommandFlag[i]._FmapdpFlag_diffusion, + numprocs,myid); + + //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed. + MPI_Barrier(MPI_COMM_WORLD); + + //Collecting data from all files generated by cores, then write to fmap_p.out + if(myid==0) + { + // system("cp 0fmap.out fmap_p.out"); + + FILE *outf; + if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL) + { + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmapdp_file); + exit_(1); + } + + FILE *fp; + char buffer[81]; + char FmapFile[50]; + for(int j=0; j<numprocs; j++) + { + FmapFile[0]='\0'; + sprintf(FmapFile,"%d",j); + strcat(FmapFile,UserCommandFlag[i]._FmapdpFlag_fmapdp_file); + + if((fp = fopen(FmapFile, "r")) == NULL) + { + fprintf(stdout, "%s: error while opening file.\n", FmapFile); + exit_(1); + } + + while(fgets(buffer, 80, fp) != NULL) + { + fputs(buffer, outf); + buffer[0]='\0'; + } + fclose(fp); + } + fclose(outf); + } + MPI_Barrier(MPI_COMM_WORLD); + printf("Process %d finished off momentum fmap.\n", myid); + + #else + fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, + UserCommandFlag[i]._FmapdpFlag_nxpoint, + UserCommandFlag[i]._FmapdpFlag_nepoint, + UserCommandFlag[i]._FmapdpFlag_nturn, + UserCommandFlag[i]._FmapdpFlag_xmax, + UserCommandFlag[i]._FmapdpFlag_emax, + UserCommandFlag[i]._FmapdpFlag_z, + UserCommandFlag[i]._FmapdpFlag_diffusion); + #endif } // // if (CodeComparaisonFlag) { @@ -539,9 +622,23 @@ int main(int argc, char *argv[]) { exit_(1); }; + /* calculate momentum acceptance*/ + printf("\n Calculate momentum acceptance: \n"); + + #if MPI_EXEC + /* calculate momentum acceptance*/ + //initialization for parallel computing + int myid=0, numprocs=0; + int namelen=0; + char processor_name[MPI_MAX_PROCESSOR_NAME]=""; + MPI_Comm_size(MPI_COMM_WORLD, &numprocs); + MPI_Comm_rank(MPI_COMM_WORLD, &myid); + MPI_Get_processor_name(processor_name, &namelen); + printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); + printf("\n Calculate momentum acceptance: \n"); - MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file, + MomentumAcceptance_p(UserCommandFlag[i]._MomentumAccFlag_momacc_file, UserCommandFlag[i]._MomentumAccFlag_istart, UserCommandFlag[i]._MomentumAccFlag_istop, UserCommandFlag[i]._MomentumAccFlag_deltaminp, @@ -551,7 +648,127 @@ int main(int argc, char *argv[]) { UserCommandFlag[i]._MomentumAccFlag_deltamaxn, UserCommandFlag[i]._MomentumAccFlag_nstepn, UserCommandFlag[i]._MomentumAccFlag_nturn, - UserCommandFlag[i]._MomentumAccFlag_zmax); + UserCommandFlag[i]._MomentumAccFlag_zmax, + numprocs,myid); + //merge the files + //Synchronize all cores, till finish cal of different region, + //then files from the cores will be processed. + MPI_Barrier(MPI_COMM_WORLD); + + //Collect data from all files generated by different cores, then write to user defined file + if(myid==0) + { + //merge the phase.out from different cpus + FILE *outf1; //phase.out, for debugging + FILE *fp1; + char buffer1[81]; + char PhaseFile[50]; + if ((outf1 = fopen("phase_p.out", "w")) == NULL) + { + fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out"); + exit_(1); + } + + for(int j=0; j<numprocs; j++) + { + PhaseFile[0]='\0'; + sprintf(PhaseFile,"%d",j); + strcat(PhaseFile,"phase.out"); + + if((fp1 = fopen(PhaseFile, "r")) == NULL){ + fprintf(stdout, "%s: error while opening file.\n", PhaseFile); + exit_(1); + } + + while(fgets(buffer1, 80, fp1) != NULL) + { + fputs(buffer1, outf1); + buffer1[0]='\0'; + } + fclose(fp1); + } + + //collect data for the momentum acceptance + FILE *outf2; //momentum acceptance + FILE *fp2; + char buffer2[81]; + char MAFile[50]; + if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL) + { + fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file); + exit_(1); + } + + //Collect data in Positive region + for(int j=0; j<numprocs; j++) + { + MAFile[0]='\0'; + sprintf(MAFile,"%d",j); + strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file); + + if((fp2 = fopen(MAFile, "r")) == NULL) + { + fprintf(stdout, "%s: error while opening file.\n", MAFile); + exit_(1); + } + + while(fgets(buffer2, 80, fp2) != NULL) + { + if(strstr(buffer2,"Negative")!=0) break; + + fputs(buffer2, outf2); + buffer2[0]='\0'; //reset buffer to NULL + } + buffer2[0]='\0'; + fclose(fp2); + } + + fprintf(outf2,"\n"); // A void line + + //Collect data in Nagative region + for(int j=0; j<numprocs; j++) + { + MAFile[0]='\0'; + sprintf(MAFile,"%d",j); + strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file); + + if((fp2 = fopen(MAFile, "r")) == NULL) + { + fprintf(stdout, "%s: error while opening file.\n", MAFile); + exit_(1); + } + + bool found=false; + while(fgets(buffer2, 80, fp2) != NULL) + { + if( (strstr(buffer2,"Negative")==0) && (!found)) { buffer2[0]='\0'; continue;} + if( (strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;} + + fputs(buffer2, outf2); + buffer2[0]='\0'; + } + buffer2[0]='\0'; + fclose(fp2); + } + fclose(outf2); + + } + MPI_Barrier(MPI_COMM_WORLD); + printf("process %d finished momentum acceptance.\n", myid); + +#else + MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file, + UserCommandFlag[i]._MomentumAccFlag_istart, + UserCommandFlag[i]._MomentumAccFlag_istop, + UserCommandFlag[i]._MomentumAccFlag_deltaminp, + UserCommandFlag[i]._MomentumAccFlag_deltamaxp, + UserCommandFlag[i]._MomentumAccFlag_nstepp, + UserCommandFlag[i]._MomentumAccFlag_deltaminn, + UserCommandFlag[i]._MomentumAccFlag_deltamaxn, + UserCommandFlag[i]._MomentumAccFlag_nstepn, + UserCommandFlag[i]._MomentumAccFlag_nturn, + UserCommandFlag[i]._MomentumAccFlag_zmax); + #endif /* restore the initial values*/ globval.Cavity_on = cavityflag; @@ -765,11 +982,16 @@ Based on parts of functions get_param( ) & ID_corr(), etc in nsls-ii_lib.cc printlatt("linlat1.out"); } - else - printf("Wrong!!!!!"); + else{ + printf("Command string %s is invalid!!!!!\n ",CommandStr); + exit_(1); + }//give an error message }//end of looking for user defined flag - + + #if MPI_EXEC + MPI_Finalize(); //finish parallel computation + #endif return 0; }//end of main()