Enable the CPU gradient flow implementation to write out the gauge fi…

…eld at specified intervals

P_M_eta.c

@@ -37,13 +37,13 @@
 #include "X_psi.h"
 #include "gamma.h"
 
-double rnorm=-1;
+double r_norm=-1;
 
-/* |R>=rnorm^2 Q^2 |S> */
+/* |R>=r_norm^2 Q^2 |S> */
 void norm_X_sqr_psi(spinor * const R, spinor * const S,
                     double const mstar);
 
-/* |R>=rnorm Q|S> */
+/* |R>=r_norm Q|S> */
 void norm_X_n_psi(spinor * const R, spinor * const S,
                   const int n, double const mstar);
 
@@ -214,7 +214,7 @@ void norm_X_sqr_psi(spinor * const R, spinor * const S, double const mstar) {
     printf("using X_psiSquare.\n");
     X_psiSquare(R, S, mstar);
   }
-  mul_r(R, rnorm*rnorm, R, VOLUME);
+  mul_r(R, r_norm*r_norm, R, VOLUME);
 
 
   free(aux_);
@@ -241,7 +241,7 @@ void norm_X_n_psi(spinor * const R, spinor * const S,
   /* Here is where we have to include our operator which in this case is
      X = 1 - (2M^2)/(D_m^dagger*D_m + M^2)  */
     X_psi(R, aux, mstar);
-    npar *= rnorm;
+    npar *= r_norm;
   }
   mul_r(R, npar, R, VOLUME);
 
@@ -312,7 +312,7 @@ void X_over_sqrt_X_sqr(spinor * const R, double * const c,
 	assign(R, aux, VOLUME);//=0
       }
       else{
-      /*|R>=rnorm^2 X^2|aux> -> since aux=d -> |R>=rnorm^2 Q^2|d>*/
+      /*|R>=r_norm^2 X^2|aux> -> since aux=d -> |R>=r_norm^2 Q^2|d>*/
         norm_X_sqr_psi(R, aux, mstar);//WARNING: - maybe we have to pass this point only when j=n-2, because R is not manipulated  in the loop body.
                                       //         - seems to setup d_n-1=0
       }
@@ -333,7 +333,7 @@ void X_over_sqrt_X_sqr(spinor * const R, double * const c,
     if(0) assign_sub_lowest_eigenvalues(R, d, no_eigenvalues-1, VOLUME);
     else assign(R, d, VOLUME);
 
-    /*|aux>=rnorm^2 Q^2|R> */
+    /*|aux>=r_norm^2 Q^2|R> */
     norm_X_sqr_psi(aux, R, mstar);
     temp1=-1.0;
     temp2=c[0]/2.;

benchmark.c

@@ -452,7 +452,7 @@ int main(int argc,char *argv[])
     printf("# Performing parallel IO test ...\n");
   }
   xlfInfo = construct_paramsXlfInfo(0.5, 0);
-  write_gauge_field( "conf.test", 64, xlfInfo);
+  write_gauge_field( "conf.test", 64, xlfInfo, g_gauge_field);
   free(xlfInfo);
   if(g_proc_id==0) {
     printf("# done ...\n");

hmc_tm.c

@@ -433,7 +433,7 @@ int main(int argc,char *argv[]) {
             fprintf(stdout, "# Writing gauge field to %s.\n", tmp_filename);
 
           xlfInfo = construct_paramsXlfInfo(plaquette_energy/(6.*VOLUME*g_nproc), trajectory_counter);
-          status = write_gauge_field( tmp_filename, gauge_precision_write_flag, xlfInfo);
+          status = write_gauge_field( tmp_filename, gauge_precision_write_flag, xlfInfo, g_gauge_field);
           free(xlfInfo);
 
           if (status) {

hopping_test.c

@@ -249,7 +249,7 @@ int main(int argc,char *argv[])
 
   random_gauge_field(reproduce_randomnumber_flag, g_gauge_field);
   if ( startoption == 2 ) {  /* restart */ 
-    write_gauge_field(gauge_input_filename,gauge_precision_write_flag,xlfInfo);
+    write_gauge_field(gauge_input_filename,gauge_precision_write_flag,xlfInfo,g_gauge_field);
   } else if ( startoption == 0 ) { /* cold */
     unit_g_gauge_field();
   } else if (startoption == 3 ) { /* continue */

io/gauge.h

@@ -33,8 +33,8 @@
 int read_gauge_field(char *filename, su3 ** const gf);
 int read_binary_gauge_data(READER *reader, DML_Checksum *checksum, paramsIldgFormat * ildgformat, su3 ** const gf);
 
-int write_gauge_field(char * filename, int prec, paramsXlfInfo const *xlfInfo);
-int write_binary_gauge_data(WRITER * writer, const int prec, DML_Checksum * checksum);
+int write_gauge_field(char * filename, int prec, paramsXlfInfo const *xlfInfo, su3 ** const gf);
+int write_binary_gauge_data(WRITER * writer, const int prec, DML_Checksum * checksum, su3 ** const gf);
 
 void write_ildg_format(WRITER *writer, paramsIldgFormat const *format);
 

io/gauge_write.c

@@ -19,7 +19,7 @@
 
 #include "gauge.ih"
 
-int write_gauge_field(char * filename, const int prec, paramsXlfInfo const *xlfInfo)
+int write_gauge_field(char * filename, const int prec, paramsXlfInfo const *xlfInfo, su3 ** const gf)
 {
   WRITER * writer = NULL;
   uint64_t bytes;
@@ -40,7 +40,7 @@ int write_gauge_field(char * filename, const int prec, paramsXlfInfo const *xlfI
 
   /* Both begin and end bit are 0, the message is begun with the format, and will end with the checksum */
   write_header(writer, 0, 0, "ildg-binary-data", bytes);
-  status = write_binary_gauge_data(writer, prec, &checksum);
+  status = write_binary_gauge_data(writer, prec, &checksum, gf);
   write_checksum(writer, &checksum, NULL);
 
   if (g_cart_id == 0 && g_debug_level > 0)

io/gauge_write_binary.c

@@ -23,7 +23,7 @@
          Probably should be done better in the future. AD. */
 
 #ifdef HAVE_LIBLEMON
-int write_binary_gauge_data(LemonWriter * lemonwriter, const int prec, DML_Checksum * checksum)
+int write_binary_gauge_data(LemonWriter * lemonwriter, const int prec, DML_Checksum * checksum, su3 ** const gf)
 {
   int x, xG, y, yG, z, zG, t, tG, status = 0;
   su3 tmp3[4];
@@ -60,10 +60,10 @@ int write_binary_gauge_data(LemonWriter * lemonwriter, const int prec, DML_Check
       for(y = 0; y < LY; y++) {
         for(x = 0; x < LX; x++) {
           rank = (DML_SiteRank) ((((tG + t)*L + zG + z)*L + yG + y)*L + xG + x);
-          memcpy(&tmp3[0], &g_gauge_field[ g_ipt[t][x][y][z] ][1], sizeof(su3));
-          memcpy(&tmp3[1], &g_gauge_field[ g_ipt[t][x][y][z] ][2], sizeof(su3));
-          memcpy(&tmp3[2], &g_gauge_field[ g_ipt[t][x][y][z] ][3], sizeof(su3));
-          memcpy(&tmp3[3], &g_gauge_field[ g_ipt[t][x][y][z] ][0], sizeof(su3));
+          memcpy(&tmp3[0], &gf[ g_ipt[t][x][y][z] ][1], sizeof(su3));
+          memcpy(&tmp3[1], &gf[ g_ipt[t][x][y][z] ][2], sizeof(su3));
+          memcpy(&tmp3[2], &gf[ g_ipt[t][x][y][z] ][3], sizeof(su3));
+          memcpy(&tmp3[3], &gf[ g_ipt[t][x][y][z] ][0], sizeof(su3));
           if(prec == 32)
             be_to_cpu_assign_double2single(filebuffer + bufoffset, tmp3, 4*sizeof(su3)/8);
           else
@@ -121,7 +121,7 @@ int write_binary_gauge_data(LemonWriter * lemonwriter, const int prec, DML_Check
 
 #else /* HAVE_LIBLEMON */
 
-int write_binary_gauge_data(LimeWriter * limewriter, const int prec, DML_Checksum * checksum)
+int write_binary_gauge_data(LimeWriter * limewriter, const int prec, DML_Checksum * checksum, su3 ** const gf)
 {
   int x, X, y, Y, z, Z, tt, t0, tag=0, id=0, status=0;
   int latticeSize[] = {T_global, g_nproc_x*LX, g_nproc_y*LY, g_nproc_z*LZ};
@@ -168,10 +168,10 @@ int write_binary_gauge_data(LimeWriter * limewriter, const int prec, DML_Checksu
             /* Rank should be computed by proc 0 only */
             rank = (DML_SiteRank) (((t0*LZ*g_nproc_z + z)*LY*g_nproc_y + y)*LX*g_nproc_x + x);
             if(g_cart_id == id) {
-              memcpy(&tmp3[0], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][1], sizeof(su3));
-              memcpy(&tmp3[1], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][2], sizeof(su3));
-              memcpy(&tmp3[2], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][3], sizeof(su3));
-              memcpy(&tmp3[3], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][0], sizeof(su3));
+              memcpy(&tmp3[0], &gf[ g_ipt[tt][X][Y][Z] ][1], sizeof(su3));
+              memcpy(&tmp3[1], &gf[ g_ipt[tt][X][Y][Z] ][2], sizeof(su3));
+              memcpy(&tmp3[2], &gf[ g_ipt[tt][X][Y][Z] ][3], sizeof(su3));
+              memcpy(&tmp3[3], &gf[ g_ipt[tt][X][Y][Z] ][0], sizeof(su3));
 
               if(prec == 32) {
                 be_to_cpu_assign_double2single(tmp2, tmp3, 4*sizeof(su3)/8);
@@ -212,10 +212,10 @@ int write_binary_gauge_data(LimeWriter * limewriter, const int prec, DML_Checksu
 #ifdef TM_USE_MPI
           else {
             if(g_cart_id == id){
-              memcpy(&tmp3[0], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][1], sizeof(su3));
-              memcpy(&tmp3[1], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][2], sizeof(su3));
-              memcpy(&tmp3[2], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][3], sizeof(su3));
-              memcpy(&tmp3[3], &g_gauge_field[ g_ipt[tt][X][Y][Z] ][0], sizeof(su3));
+              memcpy(&tmp3[0], &gf[ g_ipt[tt][X][Y][Z] ][1], sizeof(su3));
+              memcpy(&tmp3[1], &gf[ g_ipt[tt][X][Y][Z] ][2], sizeof(su3));
+              memcpy(&tmp3[2], &gf[ g_ipt[tt][X][Y][Z] ][3], sizeof(su3));
+              memcpy(&tmp3[3], &gf[ g_ipt[tt][X][Y][Z] ][0], sizeof(su3));
               if(prec == 32) {
                 be_to_cpu_assign_double2single(tmp2, tmp3, 4*sizeof(su3)/8);
                 MPI_Send((void*) tmp2, 4*sizeof(su3)/8, MPI_FLOAT, 0, tag, g_cart_grid);

meas/gradient_flow.c

@@ -33,6 +33,7 @@
 #include <string.h>
 #include <stdio.h>
 #include <math.h>
+#include <float.h>
 #include "global.h"
 #include "fatal_error.h"
 #include "aligned_malloc.h"
@@ -47,6 +48,7 @@
 #include "measure_clover_field_strength_observables.h"
 #include "meas/field_strength_types.h"
 #include "meas/measurements.h"
+#include "io/gauge.h"
 
 
 void step_gradient_flow(su3 ** x0, su3 ** x1, su3 ** x2, su3 ** z, const unsigned int type, const double eps ) {
@@ -138,7 +140,7 @@ void gradient_flow_measurement(const int traj, const int id, const int ieo) {
   measurement *meas=&measurement_list[id];
   double eps = meas->gf_eps;
   double tmax = meas->gf_tmax;
-
+  double last_write = 0.0;
 
   if( g_proc_id == 0 ) {
     printf("# Doing gradient flow measurement. id=%d\n", id);
@@ -160,10 +162,20 @@ void gradient_flow_measurement(const int traj, const int id, const int ieo) {
     fprintf(outfile, "traj t P Eplaq Esym tsqEplaq tsqEsym Wsym Qsym\n");
   }
 
-
+  // if we choose to write out the gauge field, we need extra characters 
+  // to append the trajectory id and the flow time  
+  const int gauge_config_output_filename_max_length = sizeof(meas->output_filename)/sizeof(char) +
+                                                      strlen(".000000_t0.000000");
+  char gauge_config_output_filename[gauge_config_output_filename_max_length];
+  if( meas->gf_write_interval > FLT_EPSILON && strlen(meas->output_filename) == 0 ){
+    strcpy(meas->output_filename, "flowed_conf");
+  }
 
   if (meas->external_library == QUDA_LIB) {
 #ifdef TM_USE_QUDA
+    if( meas->gf_write_interval > FLT_EPSILON ){
+      tm_debug_printf(0, 0, "WARNING: QUDA-based gradient flow, writing out flowed gauge field not yet supported!\n");
+    }
     compute_WFlow_quda( eps , tmax, traj, outfile);
 #else
     fatal_error(
@@ -199,6 +211,15 @@ void gradient_flow_measurement(const int traj, const int id, const int ieo) {
         step_gradient_flow(vt.field, x1.field, x2.field, z.field, 0, eps);
         measure_clover_field_strength_observables(vt.field, &fso[step]);
         P[step] = measure_plaquette(vt.field) / (6.0 * VOLUME * g_nproc);
+
+        if( t[step] - last_write > meas->gf_write_interval ){
+          last_write = t[step];
+          snprintf(gauge_config_output_filename, gauge_config_output_filename_max_length,
+                   "%s.%06d_t%1.6lf", meas->output_filename, traj, t[step]);
+          paramsXlfInfo *xlfInfo = construct_paramsXlfInfo(P[step], traj); 
+          write_gauge_field(gauge_config_output_filename, 64, xlfInfo, vt.field);
+          free(xlfInfo);
+        }
       }
       W = t[1] * t[1] * (2 * fso[1].E + t[1] * ((fso[2].E - fso[0].E) / (2 * eps)));
       tsqE = t[1] * t[1] * fso[1].E;

meas/measurements.h

@@ -65,9 +65,14 @@ typedef struct {
   double gf_eps;
   // maximum flow time for gf measuremnt
   double gf_tmax;
+  // interval at which the gradient flow routine should write the flowed gauge field to disk
+  double gf_write_interval;
 
   /* functions for the measurement */
   void (*measurefunc) (const int traj, const int id, const int ieo);
+
+  /* output filename, optionally used by the gradient flow measurement, for example */
+  char output_filename[500];
 
   ExternalLibrary external_library;
 } measurement;

operator/Dov_psi.c

@@ -58,20 +58,20 @@
 #include "solver/dirac_operator_eigenvectors.h"
 
 void addproj_q_invsqrt(spinor * const Q, spinor * const P, const int n, const int N);
-/* |R>=rnorm^2 Q^2 |S> */
+/* |R>=r_norm^2 Q^2 |S> */
 void norm_Q_sqr_psi(spinor * const R, spinor * const S, 
-		    const double rnorm);
-/* void norm_Q_n_psi(spinor *R, spinor *S, double m, int n, double rnorm)  */
+		    const double r_norm);
+/* void norm_Q_n_psi(spinor *R, spinor *S, double m, int n, double r_norm)  */
 /*  norm_Q_n_psi makes multiplication of any power of                      */
 /*  Q== gamma5*D_W, initial vector S, final R, finally the                 */
-/* vector R is multiplied by a factor rnorm^n                              */
-/* |R>=rnorm^n Q^n |S>  where m is a mass                                  */
+/* vector R is multiplied by a factor r_norm^n                              */
+/* |R>=r_norm^n Q^n |S>  where m is a mass                                  */
 void norm_Q_n_psi(spinor * const R, spinor * const S, 
-		  const int n, const double rnorm);
+		  const int n, const double r_norm);
 /* this is Q/sqrt(Q^2) */
 void Q_over_sqrt_Q_sqr(spinor * const R, double * const c, 
 		       const int n, spinor * const S,
-		       const double rnorm, const double minev);
+		       const double r_norm, const double minev);
 
 double ov_s = 0.6;
 double m_ov = 0.;
@@ -287,9 +287,9 @@ void addproj_q_invsqrt(spinor * const Q, spinor * const P, const int n, const in
 }
 
 
-/* |R>=rnorm^2 Q^2 |S> */
+/* |R>=r_norm^2 Q^2 |S> */
 void norm_Q_sqr_psi(spinor * const R, spinor * const S, 
-		    const double rnorm) { 
+		    const double r_norm) { 
 
   spinor *aux;
   aux=lock_Dov_WS_spinor(1);
@@ -301,19 +301,19 @@ void norm_Q_sqr_psi(spinor * const R, spinor * const S,
   gamma5(aux, R, VOLUME);
   D_psi(R, aux);
   gamma5(R, R, VOLUME);
-  mul_r(R, rnorm*rnorm, R, VOLUME);
+  mul_r(R, r_norm*r_norm, R, VOLUME);
 
   unlock_Dov_WS_spinor(1);
   return;
 }
 
-/* void norm_Q_n_psi(spinor *R, spinor *S, double m, int n, double rnorm)  */
+/* void norm_Q_n_psi(spinor *R, spinor *S, double m, int n, double r_norm)  */
 /*  norm_Q_n_psi makes multiplication of any power of                      */
 /*  Q== gamma5*D_W, initial vector S, final R, finally the                 */
-/* vector R is multiplied by a factor rnorm^n                              */
-/* |R>=rnorm^n Q^n |S>                                                     */
+/* vector R is multiplied by a factor r_norm^n                              */
+/* |R>=r_norm^n Q^n |S>                                                     */
 void norm_Q_n_psi(spinor * const R, spinor * const S, 
-		  const int n, const double rnorm) { 
+		  const int n, const double r_norm) { 
 
   int i;
   double npar = 1.;
@@ -328,7 +328,7 @@ void norm_Q_n_psi(spinor * const R, spinor * const S,
     D_psi(R, aux);
     /* Term -1-s is done in D_psi! does this comment make sense for HMC? */
     gamma5(aux, R, VOLUME);
-    npar *= rnorm;
+    npar *= r_norm;
   }
   mul_r(R, npar, aux, VOLUME);
   unlock_Dov_WS_spinor(1);
@@ -337,7 +337,7 @@ void norm_Q_n_psi(spinor * const R, spinor * const S,
 
 void Q_over_sqrt_Q_sqr(spinor * const R, double * const c, 
 		       const int n, spinor * const S,
-		       const double rnorm, const double minev) {
+		       const double r_norm, const double minev) {
 
   int j;
   double fact1, fact2, temp1, temp2, temp3, temp4, maxev, tnorm;
@@ -380,7 +380,7 @@ void Q_over_sqrt_Q_sqr(spinor * const R, double * const c,
 	assign(aux, d, VOLUME);
       }
 
-      norm_Q_sqr_psi(R, aux, rnorm);
+      norm_Q_sqr_psi(R, aux, r_norm);
       temp1=-1.0;
       temp2=c[j];
       assign_mul_add_mul_add_mul_add_mul_r(d, R, dd, aux3, fact2, fact1, temp1, temp2, VOLUME);
@@ -390,13 +390,13 @@ void Q_over_sqrt_Q_sqr(spinor * const R, double * const c,
     if(1) assign_sub_lowest_eigenvalues(R, d, no_eigenvalues-1, VOLUME);
     else assign(R, d, VOLUME);
 
-    norm_Q_sqr_psi(aux, R, rnorm);
+    norm_Q_sqr_psi(aux, R, r_norm);
     temp1=-1.0;
     temp2=c[0]/2.;
     temp3=fact1/2.;
     temp4=fact2/2.;
     assign_mul_add_mul_add_mul_add_mul_r(aux, d, dd, aux3, temp3, temp4, temp1, temp2, VOLUME);
-    norm_Q_n_psi(R, aux, 1, rnorm);
+    norm_Q_n_psi(R, aux, 1, r_norm);
   }
   else {
     /* Use the adaptive precision version using the forward recursion 
@@ -406,7 +406,7 @@ void Q_over_sqrt_Q_sqr(spinor * const R, double * const c,
     /* d = T_0(Q^2) */
     assign(d, aux3, VOLUME);
     /* dd = T_1(Q^2) */
-    norm_Q_sqr_psi(dd, d, rnorm);
+    norm_Q_sqr_psi(dd, d, r_norm);
     temp3 = fact1/2.;
     temp4 = fact2/2.;  
     assign_mul_add_mul_r(dd, d, temp3, temp4, VOLUME);
@@ -418,7 +418,7 @@ void Q_over_sqrt_Q_sqr(spinor * const R, double * const c,
     temp1=-1.0;
     for (j = 2; j <= n-1; j++) {
       /* aux = T_j(Q^2) = 2 Q^2 T_{j-1}(Q^2) - T_{j-2}(Q^2) */
-      norm_Q_sqr_psi(aux, dd, rnorm);
+      norm_Q_sqr_psi(aux, dd, r_norm);
       assign_mul_add_mul_add_mul_r(aux, dd, d, fact1, fact2, temp1, VOLUME);
       /* r = r + c_j T_j(Q^2) */
       temp2 = c[j];
@@ -444,7 +444,7 @@ void Q_over_sqrt_Q_sqr(spinor * const R, double * const c,
 
     /* r = Q r */
     assign(aux, R, VOLUME); 
-    norm_Q_n_psi(R, aux, 1, rnorm);
+    norm_Q_n_psi(R, aux, 1, r_norm);
 
   }
   /* add in piece from projected subspace */

operator/Dov_psi.h

@@ -68,7 +68,7 @@ void Qov_sq_psi_prec(spinor * const P, spinor * const S);
 
 void Q_over_sqrt_Q_sqr(spinor * const R, double * const c, 
 		       const int n, spinor * const S,
-		       const double rnorm, const double minev);
+		       const double r_norm, const double minev);
 
 void calculateOverlapPolynomial();