Merge remote-tracking branch 'origin/quda_work' into deriv_mg_tune

quda_interface.c

@@ -99,7 +99,6 @@
 #include <string.h>
 #include <math.h>
 #include <float.h>
-#include <stdbool.h> // boolean types in C
 #include "quda_interface.h"
 #include "quda_types.h"
 #include "boundary.h"
@@ -270,7 +269,6 @@ void _setDefaultQudaParam(void){
 
   inv_param.residual_type = (QudaResidualType)(QUDA_L2_RELATIVE_RESIDUAL);
   inv_param.tol_hq = 0.1;
-  // alternative reliable does not seem to work well with twisted mass (clover) fermions
   inv_param.use_alternative_reliable = 0;
 
   // Tests show that setting reliable_delta = 1e-1 results in good time to solution and good
@@ -364,7 +362,6 @@ void _setVerbosityQuda(){
 
   // general verbosity
   setVerbosityQuda(gen_verb, "# QUDA: ", stdout);
-
 }
 
 void set_force_gauge_param( QudaGaugeParam * f_gauge_param){
@@ -2932,4 +2929,57 @@ void quda_mg_tune_params(void * spinorOut, void * spinorIn, const int max_iter){
   free(tunable_params); 
 }
 
+void  compute_WFlow_quda(const double eps, const double tmax, const int traj, FILE* outfile){
+  tm_stopwatch_push(&g_timers, __func__, "");
+
+  _initQuda();
+  _loadGaugeQuda(NO_COMPRESSION);//check here the input
+
+  QudaGaugeSmearParam wflow_params = newQudaGaugeSmearParam();
+  wflow_params.smear_type = QUDA_GAUGE_SMEAR_WILSON_FLOW; 
+  wflow_params.n_steps = (int)(tmax / eps) + 3; 
+  wflow_params.epsilon = eps;
+  wflow_params.meas_interval = 1;
+
+  int n_meas= wflow_params.n_steps / wflow_params.meas_interval + 1 ;
+  QudaGaugeObservableParam *obs_param;
+  obs_param = (QudaGaugeObservableParam*) malloc(sizeof(QudaGaugeObservableParam) * n_meas);
+  for (int i=0; i<n_meas; i++){
+    obs_param[i] = newQudaGaugeObservableParam();   
+    obs_param[i].compute_plaquette = QUDA_BOOLEAN_TRUE;
+    obs_param[i].compute_qcharge = QUDA_BOOLEAN_TRUE; 
+    obs_param[i].su_project = QUDA_BOOLEAN_TRUE; 
+  }
+
+  setVerbosityQuda(QUDA_SILENT, "# QUDA: ", stdout);
+  performWFlowQuda(&wflow_params, obs_param);
+  _setVerbosityQuda();
+
+  tm_debug_printf(0, 3, "traj t P Eplaq Esym tsqEplaq tsqEsym Wsym Qsym\n");  
+
+  for(int i=1; i< wflow_params.n_steps; i+=2){  
+    const double t1 = i*eps;
+    const double P = obs_param[i].plaquette[0];
+    const double E0 = obs_param[i-1].energy[0]; // E(t=t0)
+    const double E1 = obs_param[i].energy[0]; // E(t=t1)
+    const double E2 = obs_param[i+1].energy[0]; // E(t=t2)
+    const double W = t1*t1 * (2 * E1 + t1 * ((E2 - E0) / (2 * eps)));
+    const double Q = -obs_param[i].qcharge; // topological charge Q
+
+    tm_debug_printf(0, 3,
+      "# GRADFLOW: sym(plaq)  t=%lf 1-P(t)=%1.8lf E(t)=%2.8lf(%2.8lf) t^2E=%2.8lf(%2.8lf) "
+      "W(t)=%2.8lf Q(t)=%.8lf \n",
+      t1, 1 - P, E1, 36 * (1 - P), t1*t1*E1, t1*t1 * 36 * (1 - P), W,  Q);
+
+    if (g_proc_id == 0) {
+      fprintf(outfile, "%06d %f %2.12lf %2.12lf %2.12lf %2.12lf %2.12lf %2.12lf %.12lf \n", traj,
+              t1, P, 36 * (1 - P), E1, t1 * t1 * 36 * (1 - P), t1*t1*E1, W, Q);
+      fflush(outfile);
+    }
+
+  }
+
+  free(obs_param);
+  tm_stopwatch_pop(&g_timers, 0, 1, "TM_QUDA");
+}
 

read_input.l

@@ -234,6 +234,7 @@ static inline double fltlist_next_token(int * const list_end){
   double mstarsq;
   int no_sources_z2;
 
+
   double mixcg_innereps;
   int mixcg_maxinnersolverit;
 
@@ -739,8 +740,9 @@ static inline double fltlist_next_token(int * const list_end){
 ^DflFieldIter{EQL}                 BEGIN(DFLFIELDITER);
 ^DflPolyIter{EQL}                  BEGIN(DFLPOLYITER);
 
-^BeginDeflation		                 BEGIN(INITDEFLATION);
-^BeginDDalpha       	             BEGIN(INITMULTIGRID);
+
+^BeginDeflation		           BEGIN(INITDEFLATION);
+^BeginDDalpha       	           BEGIN(INITMULTIGRID);
 ^MixCGInnerEps{EQL}                BEGIN(MIXCGEPS);
 ^MixCGMaxIter{EQL}                 BEGIN(MIXCGIT);
 
@@ -750,6 +752,7 @@ static inline double fltlist_next_token(int * const list_end){
 ^BarrierMonomialsConverge{EQL}     BEGIN(BARRIER_MONOMIALS_CONVERGE);
 
 
+
 <INITDEFLATION>Init{SPC}* {
  if(myverbose) printf("Initialising DEFLATION line %d\n", line_of_file); 
  BEGIN(DEFLATION);
@@ -3604,7 +3607,7 @@ static inline double fltlist_next_token(int * const list_end){
   BEGIN(comment_caller);
 }
 
-<INITMONOMIAL,DETMONOMIAL,CLDETMONOMIAL,CLDETRATMONOMIAL,CLDETRATRWMONOMIAL,NDPOLYMONOMIAL,NDRATMONOMIAL,NDRATCORMONOMIAL,NDCLRATMONOMIAL,NDCLRATCORMONOMIAL,CLPOLYMONOMIAL,GAUGEMONOMIAL,INTEGRATOR,INITINTEGRATOR,INITMEASUREMENT,PIONNORMMEAS,ONLINEMEAS,ORIENTEDPLAQUETTESMEAS,GRADIENTFLOWMEAS,INITOPERATOR,TMOP,DBTMOP,OVERLAPOP,WILSONOP,CLOVEROP,DBCLOVEROP,POLYMONOMIAL,PLOOP,INITGPU,GPU,RATMONOMIAL,RATCORMONOMIAL,CLRATMONOMIAL,CLRATCORMONOMIAL,INITDEFLATION,DEFLATION,INITMULTIGRID,MULTIGRID,INITEXTERNALINVERTER,QUDAINVERTER,QPHIXINVERTER,NDDETRATMONOMIAL,NDCLDETRATMONOMIAL,TUNEQUDAMGPARAMS>{SPC}*\n   {
+<INITMONOMIAL,DETMONOMIAL,CLDETMONOMIAL,CLDETRATMONOMIAL,CLDETRATRWMONOMIAL,NDPOLYMONOMIAL,NDRATMONOMIAL,NDRATCORMONOMIAL,NDCLRATMONOMIAL,NDCLRATCORMONOMIAL,CLPOLYMONOMIAL,GAUGEMONOMIAL,INTEGRATOR,INITINTEGRATOR,INITMEASUREMENT,PIONNORMMEAS,ONLINEMEAS,ORIENTEDPLAQUETTESMEAS,GRADIENTFLOWMEAS,INITOPERATOR,TMOP,DBTMOP,OVERLAPOP,WILSONOP,CLOVEROP,DBCLOVEROP,POLYMONOMIAL,PLOOP,RATMONOMIAL,RATCORMONOMIAL,CLRATMONOMIAL,CLRATCORMONOMIAL,INITDEFLATION,DEFLATION,INITMULTIGRID,MULTIGRID,INITEXTERNALINVERTER,QUDAINVERTER,QPHIXINVERTER,NDDETRATMONOMIAL,NDCLDETRATMONOMIAL,TUNEQUDAMGPARAMS>{SPC}*\n   {
   line_of_file++;
 }
 <*>{SPC}*\n                       {
@@ -3776,6 +3779,7 @@ int read_input(char * conf_file){
   mstarsq = _default_mstarsq;
   no_sources_z2 = _default_no_sources_z2;
 
+
   mixcg_innereps = _default_mixcg_innereps;
   mixcg_maxinnersolverit = _default_mixcg_maxinnersolverit;