defect added

Author: jing2li

src/lib/CMakeLists.txt

@@ -78,6 +78,12 @@ list(
   ALLOC_SRC_C
   alloc/alloc.c)
 
+list(
+  APPEND
+  PTBC_SRC_C
+  ptbc/defect.c
+)
+
 list(
   APPEND
   SOLVER_SRC_C
@@ -401,6 +407,7 @@ list(
   ${IO_SRC_C}
   ${INIT_SRC_C}
   ${APP_CONTEXT_SRC_C}
+  ${PTBC_SRC_C}
   ${ALLOC_SRC_C}
   ${SOLVER_SRC_C}
   ${TEST_SRC_C}

src/lib/get_staples.c

@@ -29,6 +29,7 @@
 #include "start.h"
 #include "su3.h"
 #include "su3adj.h"
+#include <ptbc.h>
 
 void get_staples(su3* const staple, const int x, const int mu, const su3** in_gauge_field) {
   int iy;
@@ -41,20 +42,25 @@ void get_staples(su3* const staple, const int x, const int mu, const su3** in_ga
       w1 = &in_gauge_field[x][k];
       w2 = &in_gauge_field[g_iup[x][k]][mu];
       w3 = &in_gauge_field[g_iup[x][mu]][k];
+      double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
 
       /* st = w2 * w3^d */
       _su3_times_su3d(st, *w2, *w3);
       /* v = v + w1 * st */
-      _su3_times_su3_acc(*staple, *w1, st);
+      //_su3_times_su3_acc(*staple, *w1, st);
+      _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
 
       iy = g_idn[x][k];
       w1 = &in_gauge_field[iy][k];
       w2 = &in_gauge_field[iy][mu];
       w3 = &in_gauge_field[g_iup[iy][mu]][k];
+      double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
+
       /* st = w2 * w3 */
       _su3_times_su3(st, *w2, *w3);
       /* v = v + w1^d * st */
-      _su3d_times_su3_acc(*staple, *w1, st);
+      //_su3d_times_su3_acc(*staple, *w1, st);
+      _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac1);
     }
   }
 }
@@ -71,20 +77,24 @@ void get_spacelike_staples(su3* const staple, const int x, const int mu,
       w1 = &in_gauge_field[x][k];
       w2 = &in_gauge_field[g_iup[x][k]][mu];
       w3 = &in_gauge_field[g_iup[x][mu]][k];
+      double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
 
       /* st = w2 * w3^d */
       _su3_times_su3d(st, *w2, *w3);
       /* v = v + w1 * st */
-      _su3_times_su3_acc(*staple, *w1, st);
+      //_su3_times_su3_acc(*staple, *w1, st);
+      _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
 
       iy = g_idn[x][k];
       w1 = &in_gauge_field[iy][k];
       w2 = &in_gauge_field[iy][mu];
       w3 = &in_gauge_field[g_iup[iy][mu]][k];
+      double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
       /* st = w2 * w3 */
       _su3_times_su3(st, *w2, *w3);
       /* v = v + w1^d * st */
-      _su3d_times_su3_acc(*staple, *w1, st);
+      //_su3d_times_su3_acc(*staple, *w1, st);
+      _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac1);
     }
   }
 }
@@ -101,19 +111,23 @@ void get_timelike_staples(su3* const staple, const int x, const int mu,
     w1 = &in_gauge_field[x][k];
     w2 = &in_gauge_field[g_iup[x][k]][mu];
     w3 = &in_gauge_field[g_iup[x][mu]][k];
+    double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
 
     /* st = w2 * w3^d */
     _su3_times_su3d(st, *w2, *w3);
     /* v = v + w1 * st */
-    _su3_times_su3_acc(*staple, *w1, st);
+    //_su3_times_su3_acc(*staple, *w1, st);
+    _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
 
     iy = g_idn[x][k];
     w1 = &in_gauge_field[iy][k];
     w2 = &in_gauge_field[iy][mu];
     w3 = &in_gauge_field[g_iup[iy][mu]][k];
+    double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
     /* st = w2 * w3 */
     _su3_times_su3(st, *w2, *w3);
     /* v = v + w1^d * st */
-    _su3d_times_su3_acc(*staple, *w1, st);
+    //_su3d_times_su3_acc(*staple, *w1, st);
+    _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac1);
   }
 }

src/lib/include/ptbc.h

@@ -0,0 +1,5 @@
+#include <mpi.h>
+
+
+bool is_defect(PTBCDefect *def, int const ix, int const mu);
+double get_ptbc_coeff(int const ix, int const mu);

src/lib/measure_gauge_action.c

@@ -39,6 +39,7 @@
 #include "measure_gauge_action.h"
 #include "su3.h"
 #include "su3adj.h"
+#include <ptbc.h>
 
 double measure_plaquette(const su3 *const *const gf) {
   static double res;
@@ -74,7 +75,10 @@ double measure_plaquette(const su3 *const *const gf) {
           w = &gf[ix2][mu1];
           _su3_times_su3(pr2, *v, *w);
           _trace_su3_times_su3d(ac, pr1, pr2);
-          tr = ac + kc;
+
+          // compute local parallel tempering factor: fac_1 * fac_2 * fac_3 *  fac_4
+          double const ptbc_fac = get_ptbc_coeff(ix, mu1) * get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, mu1);
+          tr = ac * ptbc_fac + kc;
           ts = tr + ks;
           tt = ts - ks;
           ks = ts;
@@ -136,7 +140,10 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
         w = &gf[ix2][0];
         _su3_times_su3(pr2, *v, *w);
         _trace_su3_times_su3d(ac, pr1, pr2);
-        ac *= (1 + lambda);
+
+        // parallel tempering factor = fac_1 * fac_2 * fac_3 *  fac_4 computed in 2 parts
+        double const ptbc_fac = get_ptbc_coeff(ix, 0) * get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, 0);
+        ac *= (1 + lambda) * ptbc_fac;
         tr = ac + kc;
         ts = tr + ks;
         tt = ts - ks;
@@ -145,17 +152,21 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
       }
       // magnetic part
       for (int mu1 = 1; mu1 < 3; mu1++) {
+        double const ptbc_fac1 = get_ptbc_coeff(ix, mu1);
+
         ix1 = g_iup[ix][mu1];
-        for (int mu2 = mu1 + 1; mu2 < 4; mu2++) {
+        for (int mu2 = mu1 + 1; mu2 < 4; mu2++) {        
           ix2 = g_iup[ix][mu2];
+          double const ptbc_fac2 = get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, mu1);
+
           v = &gf[ix][mu1];
           w = &gf[ix1][mu2];
           _su3_times_su3(pr1, *v, *w);
           v = &gf[ix][mu2];
           w = &gf[ix2][mu1];
           _su3_times_su3(pr2, *v, *w);
           _trace_su3_times_su3d(ac, pr1, pr2);
-          ac *= (1 - lambda);
+          ac *= (1 - lambda) * ptbc_fac1 * ptbc_fac2;
           tr = ac + kc;
           ts = tr + ks;
           tt = ts - ks;

src/lib/measure_rectangles.c

@@ -45,6 +45,7 @@
 #include "measure_rectangles.h"
 #include "su3.h"
 #include "su3adj.h"
+#include <ptbc.h>
 
 double measure_rectangles(const su3 **const gf) {
   static double res;
@@ -81,6 +82,7 @@ double measure_rectangles(const su3 **const gf) {
             */
             j = g_iup[i][mu];
             k = g_iup[j][nu];
+            double const ptbc_fac0 = get_ptbc_coeff(i, mu) * get_ptbc_coeff(j, nu) * get_ptbc_coeff(k, nu);
             v = &gf[i][mu];
             w = &gf[j][nu];
             _su3_times_su3(tmp, *v, *w);
@@ -95,6 +97,7 @@ double measure_rectangles(const su3 **const gf) {
             */
             j = g_iup[i][nu];
             k = g_iup[j][nu];
+            double const ptbc_fac1 = get_ptbc_coeff(i, nu) * get_ptbc_coeff(j, nu) * get_ptbc_coeff(k, mu);
             v = &gf[i][nu];
             w = &gf[j][nu];
             _su3_times_su3(tmp, *v, *w);
@@ -105,7 +108,7 @@ double measure_rectangles(const su3 **const gf) {
             _trace_su3_times_su3d(ac, pr1, pr2);
             /* 	  printf("i mu nu: %d %d %d, ac = %e\n", i, mu, nu, ac); */
             /* Kahan summation */
-            tr = ac + kc;
+            tr = ac * ptbc_fac0 * ptbc_fac1 + kc;
             ts = tr + ks;
             tt = ts - ks;
             ks = ts;

src/lib/ptbc/defect.c

@@ -0,0 +1,62 @@
+#include <stdbool.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <mpi.h>
+#include <ptbc.h>
+#include "global.h"
+
+
+// check if a link lie in defect region
+// 1. last slice of along 2. 3d cube with 0 origin
+bool is_defect(PTBCDefect *def, int const ix, int const mu) {
+  if (ix >= VOLUME) return false;
+  
+  int const global_dim[4] = {T*g_nproc_t, LX*g_nproc_t, LY*g_nproc_x, LZ*g_nproc_z};
+
+  // check start point is at second last slice and link stays on the slice
+  if (g_coord[ix][def->along] == global_dim[def->along] - 2 && mu != def->along){
+    // xf = end point of link 
+    int xf[4] = {g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3]};
+    xf[mu]  = xf[mu] + 1;
+    int dim3[3];
+    int count = 0;
+    for (int d=0; d<4; d++) {
+      if (d != def->along) {
+        dim3[count] = d;
+        count++;
+      }
+    }
+
+    for (int d=0; d<3; d++) {
+      if (g_coord[ix][dim3[d]] >= def->Ld[d] || xf[dim3[d]] >= def->Ld[d]) {
+        return false;
+      }
+    }
+    if (g_proc_id==0) printf("defect coordinate %d %d %d %d mu = %d\n", g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3], mu);
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+// get multiplying factor of parallel tempering locally (assumne no overlapping defects!)
+double get_ptbc_coeff(int const ix, int const mu) {
+  int const inst = app()->ptbc.instance_id;
+  const PTBCInstance *instance = &(app()->ptbc.instances[inst]);
+  
+  // loop over defects
+  for (int i=0; i<instance->n_coeffs; i++) {
+    PTBCDefect *def = instance->defects[i];
+    // apply coeff if within defect
+    if (is_defect(def, ix, mu)) {
+      return instance->coefficients[i];
+    }
+  }
+
+  return 1.;
+}

src/lib/su3.h

@@ -479,6 +479,17 @@ typedef struct {
   (u).c21 = (v).c20 * (w).c01 + (v).c21 * (w).c11 + (v).c22 * (w).c21; \
   (u).c22 = (v).c20 * (w).c02 + (v).c21 * (w).c12 + (v).c22 * (w).c22;
 
+#define _real_times_su3_times_su3_acc(u, v, w, r)                       \
+  (u).c00 += (r) * ((v).c00 * (w).c00 + (v).c01 * (w).c10 + (v).c02 * (w).c20); \
+  (u).c01 += (r) * ((v).c00 * (w).c01 + (v).c01 * (w).c11 + (v).c02 * (w).c21); \
+  (u).c02 += (r) * ((v).c00 * (w).c02 + (v).c01 * (w).c12 + (v).c02 * (w).c22); \
+  (u).c10 += (r) * ((v).c10 * (w).c00 + (v).c11 * (w).c10 + (v).c12 * (w).c20); \
+  (u).c11 += (r) * ((v).c10 * (w).c01 + (v).c11 * (w).c11 + (v).c12 * (w).c21); \
+  (u).c12 += (r) * ((v).c10 * (w).c02 + (v).c11 * (w).c12 + (v).c12 * (w).c22); \
+  (u).c20 += (r) * ((v).c20 * (w).c00 + (v).c21 * (w).c10 + (v).c22 * (w).c20); \
+  (u).c21 += (r) * ((v).c20 * (w).c01 + (v).c21 * (w).c11 + (v).c22 * (w).c21); \
+  (u).c22 += (r) * ((v).c20 * (w).c02 + (v).c21 * (w).c12 + (v).c22 * (w).c22);
+
 #define _su3_times_su3_acc(u, v, w)                                     \
   (u).c00 += (v).c00 * (w).c00 + (v).c01 * (w).c10 + (v).c02 * (w).c20; \
   (u).c01 += (v).c00 * (w).c01 + (v).c01 * (w).c11 + (v).c02 * (w).c21; \