enable CPU implementation of gradient flow to write out the flowed gauge

field at specified intervals

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;

quda_interface.c

@@ -2247,7 +2247,7 @@ int invert_eo_degenerate_quda(spinor * const out,
           snprintf(outname, 200, "conf_mg_refresh_fail.%.6f.%04d", g_gauge_state.gauge_id, nstore);
           paramsXlfInfo * xlfInfo = construct_paramsXlfInfo(
               measure_plaquette((const su3**)g_gauge_field)/(6.*VOLUME*g_nproc), nstore);
-          int status = write_gauge_field(outname, 64, xlfInfo);
+          int status = write_gauge_field(outname, 64, xlfInfo, g_gauge_field);
           free(xlfInfo);
 
           char errmsg[200];

read_input.l

@@ -3063,6 +3063,8 @@ static inline double fltlist_next_token(int * const list_end){
 
   meas->gf_eps = _default_gf_eps;
   meas->gf_tmax = _default_gf_tmax;
+  meas->gf_write_interval = 0.0;
+  strcpy(meas->output_filename, "");
 
   if(!reread) {
     if(add_measurement(meas->type) < 0) {
@@ -3132,6 +3134,37 @@ static inline double fltlist_next_token(int * const list_end){
     meas->gf_tmax = c;
     if(myverbose) printf("  Maximum gradient flow time size set to %lf line %d, measurement id=%d\n", meas->gf_tmax, line_of_file, meas->id);
   }
+  {SPC}*GaugeConfigOutputFile{EQL}{FILENAME} {
+    char * argpos = strchr(yytext, '=');
+    if( argpos == NULL ){
+      char error_message[200];
+      snprintf(error_message, 200,
+              "Error in parsing of GaugeConfigOutputFile in line %d, measurement id=%d!\n",
+              line_of_file, meas->id);
+      fatal_error(error_message, __func__);
+    }
+    // we check length, allowing for some extraneous whitespace
+    if( strlen(argpos) >= (sizeof(meas->output_filename)/sizeof(char) + 20) ){
+      char error_message[200];
+      snprintf(error_message, 200,
+               "Filename passed as GaugeConfigOutputFile too long, line %d, measurement id=%d "
+               "(see read_input.l and meas/measurement.h)\n",
+               line_of_file, meas->id);
+      fatal_error(error_message, __func__); 
+    }
+    // this is not 100% safe but should be okay
+    sscanf(argpos, "= %s", meas->output_filename);
+    if(myverbose){
+      printf("  Gradient flow GaugeConfigOutputFile filename set to %s, line %d, measurement id=%d\n",
+             meas->output_filename, line_of_file, meas->id);
+    }
+  }
+  {SPC}*GaugeConfigOutputInterval{EQL}{FLT} {
+    sscanf(yytext, " %[2a-zA-Z] = %lf", name, &c);
+    meas->gf_write_interval = c;
+    if(myverbose) printf("  Gradient flow gauge config output interval set to %lf line %d, measurement id=%d\n",
+                         meas->gf_write_interval, line_of_file, meas->id);
+  }
 } 
 
 <PLOOP>{

solver/monomial_solve.c

@@ -95,7 +95,7 @@ void solve_fail_write_config_and_abort(const char * const solver) {
   char outname[200];
   snprintf(outname, 200, "conf_monomial_solve_fail.%.6f.%04d", g_gauge_state.gauge_id, nstore);
   paramsXlfInfo * xlfInfo = construct_paramsXlfInfo(measure_plaquette((const su3**)g_gauge_field)/(6.*VOLUME*g_nproc), nstore);
-  int status = write_gauge_field(outname, 64, xlfInfo);
+  int status = write_gauge_field(outname, 64, xlfInfo, g_gauge_field);
   free(xlfInfo);
   fatal_error("Error: solver reported -1 iterations.", solver);
 }

update_tm.c

@@ -193,7 +193,7 @@ int update_tm(double *plaquette_energy, double *rectangle_energy,
       if(g_proc_id == 0 && g_debug_level > 0) {
         fprintf(stdout, "# Writing gauge field to file %s.\n", tmp_filename);
       }
-      if((iostatus = write_gauge_field( tmp_filename, 64, xlfInfo) != 0 )) {
+      if((iostatus = write_gauge_field( tmp_filename, 64, xlfInfo, g_gauge_field) != 0 )) {
         /* Writing failed directly */
         fprintf(stderr, "Error %d while writing gauge field to %s\nAborting...\n", iostatus, tmp_filename);
         exit(-2);