From 04d38e296ead3bfcb6b6a05535a84d3e0a6bcfc4 Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Mon, 20 May 2024 16:31:13 +0100 Subject: [PATCH 01/17] add code for APE and Wuppertal smearing (does not compile - old complex numbers) --- Include/global.h | 17 ++ Include/observables.h | 21 +++ Include/representation.h | 2 + Include/utils.h | 12 ++ LibHR/Observables/avr_plaquette.c | 172 +++++++++++++++++ LibHR/Observables/measure_mesons.c | 170 +++++++++++++++++ LibHR/Observables/meson_measurements.c | 73 +++++++ LibHR/Observables/polyakov.c | 251 ++++++++++++++++++++++++- LibHR/Observables/sources.c | 186 ++++++++++++++++++ LibHR/Update/representation.c | 75 ++++++++ LibHR/Utils/APE_smearing.c | 128 +++++++++++++ LibHR/Utils/boundary_conditions.c | 128 +++++++++++++ LibHR/Utils/single_double_utils.c | 18 ++ LibHR/Utils/suN_utils.c | 192 +++++++++++++++++++ Spectrum/measure_spectrum.c | 33 +++- 15 files changed, 1476 insertions(+), 2 deletions(-) create mode 100644 LibHR/Utils/APE_smearing.c diff --git a/Include/global.h b/Include/global.h index 0d3b2e90b..b7da350a5 100644 --- a/Include/global.h +++ b/Include/global.h @@ -143,6 +143,15 @@ GLB_VAR(suNg_field_flt,*u_gauge_flt,=NULL); GLB_VAR(suNf_field,*u_gauge_f,=NULL); GLB_VAR(suNg_field,*u_gauge_s,=NULL); GLB_VAR(suNf_field_flt,*u_gauge_f_flt,=NULL); + +/*APE smearing*/ +GLB_VAR(suNg_field,*u_gauge_APE,=NULL); +GLB_VAR(suNg_scalar_field,*u_scalar_APE,=NULL); +GLB_VAR(suNg_field_flt,*u_gauge_APE_flt,=NULL); +GLB_VAR(suNf_field,*u_gauge_APE_f,=NULL); +GLB_VAR(suNg_field,*u_gauge_APE_s,=NULL); +GLB_VAR(suNf_field_flt,*u_gauge_APE_f_flt,=NULL); + #if defined(GAUGE_SPN) && defined(REPR_FUNDAMENTAL) GLB_VAR(suNffull_field,*cl_term,=NULL); GLB_VAR(suNffull_field,*cl_force,=NULL); @@ -161,6 +170,14 @@ GLB_VAR(int,gauge_field_active,=0); // whether gauge field interactions is activ #define pu_gauge_f(ix,mu) ((u_gauge_f->ptr)+coord_to_index(ix,mu)) #define pu_gauge_f_flt(ix,mu) ((u_gauge_f_flt->ptr)+coord_to_index(ix,mu)) +/*APE smearing*/ +#define pu_gauge_APE(ix,mu) ((u_gauge_APE->ptr)+coord_to_index(ix,mu)) +#define pu_scalar_APE(ix) ((u_scalar_APE->ptr)+ix) +#define pu_gauge_APE_flt(ix,mu) ((u_gauge_APE_flt->ptr)+coord_to_index(ix,mu)) +#define pu_gauge_APE_f(ix,mu) ((u_gauge_APE_f->ptr)+coord_to_index(ix,mu)) +#define pu_gauge_APE_f_flt(ix,mu) ((u_gauge_APE_f_flt->ptr)+coord_to_index(ix,mu)) +#define pu_gauge_tmp(ix,mu) ((u_gauge_tmp->ptr)+coord_to_index(ix,mu)) + /* input parameters */ #include "input_par.h" GLB_VAR(input_glb,glb_var,=init_input_glb(glb_var)); diff --git a/Include/observables.h b/Include/observables.h index d7a9bef87..b986607f1 100644 --- a/Include/observables.h +++ b/Include/observables.h @@ -296,4 +296,25 @@ void init_triplet_discon_correlators(); void ff_observables(); +//APE smearing +double plaq_APE(int ix,int mu,int nu); +void cplaq_APE(complex *ret,int ix,int mu,int nu); +double avr_spacial_plaquette_APE(); +void full_plaquette_APE(); +void polyakov_APE(); + +/* source smearing*/ +void smearing_function(spinor_field *source, int tau, int color, double epsilon); +void smearing_function_with_APE(spinor_field *source, int tau, int color, double epsilon); +void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear); +void create_smeared_source_with_APE(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear); + +/* Sink smearing*/ + +void smeared_propagator(spinor_field* psi, int nm , double epsilon); +void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon); + +/* Measurements*/ +void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N); + #endif diff --git a/Include/representation.h b/Include/representation.h index 328d8504b..a7b88d129 100644 --- a/Include/representation.h +++ b/Include/representation.h @@ -21,4 +21,6 @@ void _group_represent2_flt(suNf_flt* v, suNg_flt *u); void represent_gauge_field(); void represent_gauge_field_measure(); +/* APE smearing*/ +void represent_gauge_field_APE(); #endif diff --git a/Include/utils.h b/Include/utils.h index 98d3d9437..fdda4d801 100644 --- a/Include/utils.h +++ b/Include/utils.h @@ -44,6 +44,7 @@ void init_plaq_open_BCs(double *plaq_weight, double *rect_weight, double ct, dou void free_BCs(); void apply_BCs_on_represented_gauge_field(); +void apply_BCs_on_represented_gauge_field_APE(); void apply_BCs_on_fundamental_gauge_field(); void apply_BCs_on_momentum_field(suNg_av_field *force); void apply_BCs_on_spinor_field(spinor_field *sp); @@ -176,4 +177,15 @@ int reserve_wrk_space_with_pointers(suNg_field **g_wrk_out, int **i_up_wrk_out, void release_wrk_space(int id_release); void free_wrk_space(); +/* APE smearing */ +void APE_smearing(double smear_val, int Nsmear); +void assign_ud2u_APE_f(void); +#ifdef GAUGE_SPN +void cooling_SPN(suNg* g_out, suNg* g_in, suNg* staple, int cooling); +void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11); +void vmxsu2(int i1, int i2, suNgfull* A, complex B[4]); +void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11); +void vmxsu2_tau(int n1, int n2, suNgfull* A, complex B[4]); +#endif + #endif diff --git a/LibHR/Observables/avr_plaquette.c b/LibHR/Observables/avr_plaquette.c index 0a47117b1..1d8994c57 100644 --- a/LibHR/Observables/avr_plaquette.c +++ b/LibHR/Observables/avr_plaquette.c @@ -416,3 +416,175 @@ double complex avr_plaquette_wrk() return pa; } + + +double plaq_APE(int ix,int mu,int nu) +{ + int iy,iz; + double p; + suNg *v1,*v2,*v3,*v4,w1,w2,w3; + + iy=iup(ix,mu); + iz=iup(ix,nu); + + v1=pu_gauge_APE(ix,mu); + v2=pu_gauge_APE(iy,nu); + v3=pu_gauge_APE(iz,mu); + v4=pu_gauge_APE(ix,nu); + + _suNg_times_suNg(w1,(*v1),(*v2)); + _suNg_times_suNg(w2,(*v4),(*v3)); + _suNg_times_suNg_dagger(w3,w1,w2); + + _suNg_trace_re(p,w3); + +#ifdef PLAQ_WEIGHTS + if(plaq_weight==NULL) return p; + return plaq_weight[ix*16+mu*4+nu]*p; +#else + return p; +#endif +} + +void cplaq_APE(complex *ret,int ix,int mu,int nu) +{ + int iy,iz; + suNg *v1,*v2,*v3,*v4,w1,w2,w3; + + iy=iup(ix,mu); + iz=iup(ix,nu); + + v1=pu_gauge_APE(ix,mu); + v2=pu_gauge_APE(iy,nu); + v3=pu_gauge_APE(iz,mu); + v4=pu_gauge_APE(ix,nu); + + _suNg_times_suNg(w1,(*v1),(*v2)); + _suNg_times_suNg(w2,(*v4),(*v3)); + _suNg_times_suNg_dagger(w3,w1,w2); + + _suNg_trace_re(ret->re,w3); +#ifdef GAUGE_SON + ret->im=0; +#else + _suNg_trace_im(ret->im,w3); +#endif + +#ifdef PLAQ_WEIGHTS + if(plaq_weight!=NULL) { + ret->re *= plaq_weight[ix*16+mu*4+nu]; + ret->im *= plaq_weight[ix*16+mu*4+nu]; + } +#endif + +} + + +double avr_spacial_plaquette_APE() +{ + double pa=0.; + + _PIECE_FOR(&glattice,ixp) { + if(ixp==glattice.inner_master_pieces) { + _OMP_PRAGMA( master ) + /* wait for gauge field to be transfered */ + complete_gf_sendrecv(u_gauge_APE); + _OMP_PRAGMA( barrier ) + } + _SITE_FOR_SUM(&glattice,ixp,ix,pa) { + + pa+=plaq_APE(ix,2,1); + pa+=plaq_APE(ix,3,1); + pa+=plaq_APE(ix,3,2); + } + } + + global_sum(&pa, 1); + +#ifdef BC_T_OPEN + pa /= 3.0*NG*GLB_VOLUME*(GLB_T-1)/GLB_T; +#else + pa /= 3.0*NG*GLB_VOLUME; +#endif + + return pa; + +} + +void full_plaquette_APE() +{ + complex pa[6]; + double r0re = 0, r0im = 0; + double r1re = 0, r1im = 0; + double r2re = 0, r2im = 0; + double r3re = 0, r3im = 0; + double r4re = 0, r4im = 0; + double r5re = 0, r5im = 0; + + _PIECE_FOR(&glattice,ixp) + { + if(ixp == glattice.inner_master_pieces) + { + _OMP_PRAGMA( master ) + /* wait for gauge field to be transfered */ + complete_gf_sendrecv(u_gauge_APE); + _OMP_PRAGMA( barrier ) + } + + _SITE_FOR_SUM(&glattice,ixp,ix,r0re,r0im,r1re,r1im,r2re,r2im,r3re,r3im,r4re,r4im,r5re,r5im) + { + complex tmp; + + cplaq_APE(&tmp,ix,1,0); + r0re += tmp.re; + r0im += tmp.im; + + cplaq_APE(&tmp,ix,2,0); + r1re += tmp.re; + r1im += tmp.im; + + cplaq_APE(&tmp,ix,2,1); + r2re += tmp.re; + r2im += tmp.im; + + cplaq_APE(&tmp,ix,3,0); + r3re += tmp.re; + r3im += tmp.im; + + cplaq_APE(&tmp,ix,3,1); + r4re += tmp.re; + r4im += tmp.im; + + cplaq_APE(&tmp,ix,3,2); + r5re += tmp.re; + r5im += tmp.im; + + } + } + + pa[0].re=r0re; pa[0].im=r0im; + pa[1].re=r1re; pa[1].im=r1im; + pa[2].re=r2re; pa[2].im=r2im; + pa[3].re=r3re; pa[3].im=r3im; + pa[4].re=r4re; pa[4].im=r4im; + pa[5].re=r5re; pa[5].im=r5im; + + global_sum((double*)pa,12); + for(int k = 0; k < 6; k++) + { +#ifdef BC_T_OPEN + pa[k].re /= NG*GLB_VOLUME*(GLB_T-1)/GLB_T; + pa[k].im /= NG*GLB_VOLUME*(GLB_T-1)/GLB_T; +#else + pa[k].re /= NG*GLB_VOLUME; + pa[k].im /= NG*GLB_VOLUME; +#endif + } + + lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",1,0,pa[0].re,pa[0].im); + lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",2,0,pa[1].re,pa[1].im); + lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",2,1,pa[2].re,pa[2].im); + lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",3,0,pa[3].re,pa[3].im); + lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",3,1,pa[4].re,pa[4].im); + lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",3,2,pa[5].re,pa[5].im); +} \ No newline at end of file diff --git a/LibHR/Observables/measure_mesons.c b/LibHR/Observables/measure_mesons.c index 378eb0622..8cf6ee93f 100644 --- a/LibHR/Observables/measure_mesons.c +++ b/LibHR/Observables/measure_mesons.c @@ -623,3 +623,173 @@ void print_mesons(meson_observable *mo, double norm, int conf, int nm, double *m print_corr(mo, lt, conf, nm, mass, label, n_mom); zero_corrs(mo); } + + +/* Wuppertal smearing*/ + +void smeared_propagator(spinor_field* psi, int nm , double epsilon){ + int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; + double norm_factor = 1./(1.+6.*epsilon); + suNf_propagator sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down; + _propagator_zero(sp_smeared); + + spinor_field* smeared_psi = alloc_spinor_field_f(4*nm*NF,&glattice); + + for(i=0; ic[alpha].c[b].re = sp_smeared.c[b].c[alpha].c[beta].c[a].re; + _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b].im = sp_smeared.c[b].c[alpha].c[beta].c[a].im; + + } + } + + } + } + } + for (int i=0;i<4*nm*NF;i++){ + spinor_field_copy_f(psi + i, smeared_psi + i); + } + for (int i=0;i<4*nm*NF;i++){ + start_sf_sendrecv(psi + i); + complete_sf_sendrecv(psi + i); + } + free_spinor_field_f(smeared_psi); +} + +void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon){ + + int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; + double norm_factor = 1./(1.+6.*epsilon); + suNf_propagator sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down; + _propagator_zero(sp_smeared); + + spinor_field* smeared_psi = alloc_spinor_field_f(4*nm*NF, &glattice); + + for(i=0; ic[alpha].c[b].re = sp_smeared.c[b].c[alpha].c[beta].c[a].re; + _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b].im = sp_smeared.c[b].c[alpha].c[beta].c[a].im; + + } + } + } + } + } + for (int i=0;i<4*nm*NF;i++){ + spinor_field_copy_f(psi + i, smeared_psi + i); + } + for (int i=0;i<4*nm*NF;i++){ + start_sf_sendrecv(psi + i); + complete_sf_sendrecv(psi + i); + } + free_spinor_field_f(smeared_psi); +} + diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index 52fe5e328..ad23e21d7 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -761,3 +761,76 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m free_gfield_f(u_gauge_old); free_propagator_eo(); } + + + +/* Updates of Wuppertal smearing source in 2020*/ + +void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N){ + + char label[20]; + spinor_field* source = alloc_spinor_field_f(4,&glattice); + spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice); + + init_propagator_eo(nm,m,precision); + + if (APE_N != 0){ + APE_smearing(APE_epsilon, APE_N); + represent_gauge_field_APE(); + } + + int k; + for (k=0;k1){ + measure_point_mesons_momenta(meson_correlators,prop+4*k, source, nm, t, n_mom); + } + else{ + measure_mesons(meson_correlators,prop+4*k, source, nm, t); + } + } + + sprintf(label, "source_N%d_sink_N%d",Nsmear_source, 0); + print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); + + int N_diff = 10; // measure sink smearing every 10 steps + for (int N = 0; N1){ + measure_point_mesons_momenta(meson_correlators,prop+4*c, source, nm, t, n_mom); + } + else{ + measure_mesons(meson_correlators,prop+4*c, source, nm, t); + } + } + sprintf(label, "source_N%d_sink_N%d",Nsmear_source, N+N_diff); + print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); + } + + free_propagator_eo(); + free_spinor_field_f(source); + free_spinor_field_f(prop); +} \ No newline at end of file diff --git a/LibHR/Observables/polyakov.c b/LibHR/Observables/polyakov.c index 7c72c93e1..ae34d4075 100644 --- a/LibHR/Observables/polyakov.c +++ b/LibHR/Observables/polyakov.c @@ -458,4 +458,253 @@ void polyakov() #endif #undef MPIRET -#undef _MIN \ No newline at end of file +#undef _MIN + + +void polyakov_APE() { + int x[4], i4d, i3d, size3d; + int loc[4]={T,X,Y,Z}; + suNg *p; + suNg *lp; + suNg *bp; + suNg tmp; + complex poly; + double adjpoly; + double dtmp; +#ifdef WITH_MPI + int np[4]={NP_T,NP_X,NP_Y,NP_Z}; + int mpiret; +#endif /* WITH_MPI */ + + for(int mu=0; mu<4; mu++){ + + size3d=T*X*Y*Z/loc[mu]; + p=malloc(sizeof(suNg)*size3d); + lp=malloc(sizeof(suNg)*size3d); + bp=malloc(sizeof(suNg)*size3d); + + + /* local wilson lines */ + + i3d=0; + for(x[(mu+1)%4]=0; x[(mu+1)%4]c[spin] = spinor_smeared.c[spin]; + + } + } + + for (int beta=0;beta<4;++beta){ + spinor_field_copy_f(source + beta, smeared_source +beta); + } + } + + for (int beta=0;beta<4;++beta){ + start_sf_sendrecv(source + beta); + complete_sf_sendrecv(source + beta); + } + free_spinor_field_f(smeared_source); +} + +void smearing_function_with_APE(spinor_field *source, int tau, int color, double epsilon){ + + int ix, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; + double norm_factor = 1./(1.+6.*epsilon); + suNf_spinor spinor_OG, spinor_smeared, spinor_tmp, spinor_right, spinor_left, spinor_front, spinor_back, spinor_up, spinor_down; + _spinor_zero_f(spinor_OG);_spinor_zero_f(spinor_tmp);_spinor_zero_f(spinor_right);_spinor_zero_f(spinor_left); + _spinor_zero_f(spinor_front);_spinor_zero_f(spinor_back);_spinor_zero_f(spinor_up);_spinor_zero_f(spinor_down); + + spinor_field* smeared_source = alloc_spinor_field_f(4, &glattice); + for (int beta=0;beta<4;++beta){ + spinor_field_zero_f(&smeared_source[beta]); + } + + represent_gauge_field_APE(); + if(COORD[0]==tau/T){ + for (x=0; x +#include +#include +#include +#include "global.h" +#include "utils.h" +#include "suN.h" +#include "memory.h" +#include "communications.h" +#include "logger.h" +#include "complex.h" +#include "random.h" +#include "representation.h" +#include "observables.h" +#include "update.h" + +void APE_smearing(double smear_val, int Nsmear){ + + suNg staple, tr1, tr2; + int ix, iy, iz, it; + int mu, nu, mid, midpmu, midpnu, midmnu, midpmumnu; + suNg v, vout, vtmp; + + double sm1 = 1. - smear_val; + double sm2 = smear_val / 6. ; + + suNg_field *u_gauge_tmp = alloc_gfield(&glattice); + u_gauge_APE = alloc_gfield(&glattice); + + #ifdef ALLOCATE_REPR_GAUGE_FIELD + u_gauge_APE_f=alloc_gfield_f(&glattice); + #endif + + for (it=0;it_s = %1.6f\n",avr_spacial_plaquette_APE()); + + for (int N=0; N_s = %1.6f\n", Nsmear, avr_spacial_plaquette_APE()); + polyakov_APE(); + full_plaquette_APE(); + lprintf("APE",0,"APE smearing END \n"); + free_gfield(u_gauge_tmp); +} \ No newline at end of file diff --git a/LibHR/Utils/boundary_conditions.c b/LibHR/Utils/boundary_conditions.c index 15869eb60..4e11daf68 100644 --- a/LibHR/Utils/boundary_conditions.c +++ b/LibHR/Utils/boundary_conditions.c @@ -337,6 +337,33 @@ void apply_BCs_on_clover_term(suNfc_field *cl) } #endif //defined(GAUGE_SPN) && defined(REPR_FUNDAMENTAL) +static void sp_T_antiperiodic_BCs_APE(); +static void sp_X_antiperiodic_BCs_APE(); +static void sp_Y_antiperiodic_BCs_APE(); +static void sp_Z_antiperiodic_BCs_APE(); +/*static void sp_spatial_theta_BCs(double theta);*/ +static void chiSF_ds_BT_APE(double ds); + +void apply_BCs_on_represented_gauge_field_APE() { +#ifdef BC_T_ANTIPERIODIC + sp_T_antiperiodic_BCs_APE(); +#endif +#ifdef BC_X_ANTIPERIODIC + sp_X_antiperiodic_BCs_APE(); +#endif +#ifdef BC_Y_ANTIPERIODIC + sp_Y_antiperiodic_BCs_APE(); +#endif +#ifdef BC_Z_ANTIPERIODIC + sp_Z_antiperiodic_BCs_APE(); +#endif +#ifdef ROTATED_SF +#ifndef ALLOCATE_REPR_GAUGE_FIELD +#error The represented gauge field must be allocated!!! +#endif + chiSF_ds_BT_APE(BCs_pars.chiSF_boundary_improvement_ds); +#endif +} /***************************************************************************/ /* BOUNDARY CONDITIONS TO BE APPLIED ON THE REPRESENTED GAUGE FIELD */ @@ -816,6 +843,107 @@ static void gf_open_BCs() } } #endif +/***************************************************************************/ +/* BOUNDARY CONDITIONS TO BE APPLIED ON THE REPRESENTED APE GAUGE FIELD */ +/***************************************************************************/ + +static void sp_T_antiperiodic_BCs_APE() { + if(COORD[0]==0) { + int index; + int ix,iy,iz; + suNf *u; + for (ix=0;ixptr); + f = (float*)(u_gauge_APE_f_flt->ptr); + +_OMP_PRAGMA( _omp_parallel ) +_OMP_PRAGMA( _omp_for ) + for(int i=0; i<4*glattice.gsize_gauge*(sizeof(suNf)/sizeof(double)); i++) { + *(f+i) = (float)(*(d+i)); + } + } +} \ No newline at end of file diff --git a/LibHR/Utils/suN_utils.c b/LibHR/Utils/suN_utils.c index b7ef83e4a..d934e461b 100644 --- a/LibHR/Utils/suN_utils.c +++ b/LibHR/Utils/suN_utils.c @@ -447,3 +447,195 @@ void covariant_project_to_suNg(suNg *u) _suNg_times_suNg(*u, tmp1, tmp); } + +void cooling_SPN(suNg* g_out, suNg* g_in, suNg* g_tilde, int cooling){ + + suNgfull B, U_tilde, S; + _suNg_expand( U_tilde, *g_tilde); + _suNg_expand( S, *g_in); + + //_suNffull_unit(S); + + /* + lprintf("COOLING",0," U_tilde= \n",0); + for (int i=0;ic[i].re = S.c[i].re; + g_out->c[i].im = S.c[i].im; + } +} + +void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11){ + + complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array + double UMAG; // [ 0 1 ] + int i1,i2,i3,i4; // [ 2 3 ] --> [0, 1, 2, 3] + + i1 = i1col + i1col*NG; + i2 = i2col + i2col*NG; + i3 = i2col + i1col*NG; + i4 = i1col + i2col*NG; + + _complex_add_star(F11, B11->c[i1], B11->c[i2]); + _complex_mulr(F11, 0.5, F11); + _complex_sub_star(F12, B11->c[i3], B11->c[i4]); + _complex_mulr(F12, 0.5, F12); + + + _complex_mul_star(ztmp1, F11, F11); + _complex_mul_star(ztmp2, F12, F12); + UMAG = sqrt(ztmp1.re+ztmp2.re); + UMAG = 1./UMAG; + + _complex_mulr(F11, UMAG, F11); + _complex_mulr(F12, UMAG, F12); + + _complex_star(A11[0], F11); + _complex_star(A11[1], F12); + _complex_mulr(A11[2], -1., F12); + _complex_mulr(A11[3], 1., F11); + + vmxsu2(i1col,i2col,C11,A11); + vmxsu2(i1col,i2col,B11,A11); +} + +void vmxsu2(int i1, int i2, suNgfull* A, complex B[4]){ + + complex C[NG*2]; + complex ztmp1,ztmp2; + + for (int i=0; ic[i + i1*NG], B[0]); + _complex_mul(ztmp2, A->c[i + i2*NG], B[2]); + _complex_add(C[i], ztmp1, ztmp2); + + _complex_mul(ztmp1, A->c[i + i1*NG], B[1]); + _complex_mul(ztmp2, A->c[i + i2*NG], B[3]); + _complex_add(C[i+NG], ztmp1, ztmp2); + } + + for (int i=0; ic[i + i1*NG] = C[i]; + A->c[i + i2*NG] = C[i+NG]; + } +} + +void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11){ + + complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array + double UMAG; // [ 0 1 ] + int i1,i2,i3,i4; // [ 2 3 ] --> [0, 1, 2, 3] + + i1 = n1 + n1*NG; + i2 = (NG/2) + n2 + (NG/2 + n2)*NG; + i3 = NG/2 + n2 + n1*NG; + i4 = (NG/2 + n2)*NG + n1; + + _complex_add_star(F11, B11->c[i1], B11->c[i2]); + _complex_mulr(F11, 0.5, F11); + _complex_sub_star(F12, B11->c[i3], B11->c[i4]); + _complex_mulr(F12, 0.5, F12); + + _complex_mul_star(ztmp1, F11, F11); + _complex_mul_star(ztmp2, F12, F12); + UMAG = sqrt(ztmp1.re+ztmp2.re); + UMAG = 1./UMAG; + + _complex_mulr(F11, UMAG, F11); + _complex_mulr(F12, UMAG, F12); + + _complex_star(A11[0], F11); + _complex_star(A11[1], F12); + _complex_mulr(A11[2], -1., F12); + _complex_mulr(A11[3], 1., F11); + + vmxsu2_tau(n1, n2, C11, A11); + vmxsu2_tau(n1, n2, B11, A11); +} + +void vmxsu2_tau(int n1, int n2, suNgfull* A, complex B[4]){ + + suNgfull C; + complex ztmp1, ztmp2; + + _suNgfull_mul(C, 1., *A); + + for (int i=0; ic[i + n1*NG], B[0]); + _complex_mul(ztmp2, A->c[i +(NG/2+n2)*NG], B[2]); + _complex_add(C.c[i + n1*NG], ztmp1, ztmp2); + + _complex_mul(ztmp1, A->c[i + n2*NG], B[0]); + _complex_mul(ztmp2, A->c[i + (NG/2+n1)*NG], B[2]); + _complex_add(C.c[i + n2*NG], ztmp1, ztmp2); + + _complex_mul(ztmp1, A->c[i + n2*NG], B[1]); + _complex_mul(ztmp2, A->c[i + (NG/2+n1)*NG], B[3]); + _complex_add(C.c[i + (NG/2 + n1)*NG], ztmp1, ztmp2); + + _complex_mul(ztmp1, A->c[i + n1*NG], B[1]); + _complex_mul(ztmp2, A->c[i + (NG/2+n2)*NG], B[3]); + _complex_add(C.c[i + (NG/2 + n2)*NG], ztmp1, ztmp2); + } + _suNgfull_mul(*A, 1., C); +} diff --git a/Spectrum/measure_spectrum.c b/Spectrum/measure_spectrum.c index c8879478c..a8d5f5d2f 100755 --- a/Spectrum/measure_spectrum.c +++ b/Spectrum/measure_spectrum.c @@ -74,13 +74,24 @@ typedef struct _input_mesons double csw; double rho_s; double rho_t; + int source_t; + int source_x; + int source_y; + int source_z; + int smearing_ss; + double smear_epsilon_source; + int smear_N_source; + double smear_epsilon_sink; + int smear_N_sink; + int APE_N; + double APE_epsilon; //Currently only implemented for ff int nhits_hopping; //Multiplies the number of hits in the fast part of the hopping parameter expansion int degree_hopping; // The degree of the hopping parameter expasion /* for the reading function */ - input_record_t read[33]; + input_record_t read[44]; } input_mesons; #define init_input_mesons(varname) \ @@ -117,6 +128,17 @@ typedef struct _input_mesons {"hopping expansion degree", "mes:degree_hopping = %d", INT_T, &(varname).degree_hopping}, \ {"hopping expansion hits", "mes:nhits_hopping = %d", INT_T, &(varname).nhits_hopping}, \ {"Configuration list:", "mes:configlist = %s", STRING_T, &(varname).configlist}, \ + {"Smearing Source z", "mes:source_t = %d",INT_T, &(varname).source_t}, \ + {"Smearing Source x", "mes:source_x = %d",INT_T, &(varname).source_x}, \ + {"Smearing Source y", "mes:source_y = %d",INT_T, &(varname).source_y}, \ + {"Smearing Source z", "mes:source_z = %d",INT_T, &(varname).source_z}, \ + {"Smearing a lot at once", "mes:smearing_ss = %d",INT_T, &(varname).smearing_ss}, \ + {"smearing source step size", "mes:smear_epsilon_source = %lf", DOUBLE_T, &(varname).smear_epsilon_source}, \ + {"Number of smeared source steps", "mes:smear_N_source = %d",INT_T, &(varname).smear_N_source}, \ + {"smearing sink step size", "mes:smear_epsilon_sink = %lf", DOUBLE_T, &(varname).smear_epsilon_sink}, \ + {"Number of smeared sink steps", "mes:smear_N_sink = %d",INT_T, &(varname).smear_N_sink}, \ + {"Number of APE smearing steps", "mes:APE_N = %d",INT_T, &(varname).APE_N}, \ + {"APE smearing step size", "mes:APE_epsilon = %lf",DOUBLE_T, &(varname).APE_epsilon}, \ {NULL, NULL, INT_T, NULL} \ } \ } @@ -208,6 +230,12 @@ int main(int argc, char *argv[]) { lprintf("MAIN", 0, "Point sources\n"); } + if (mes_var.APE_N != 0){ + lprintf("MAIN",0,"APE smearing activate on Gaussian smearing\n"); + } + if (mes_var.smearing_ss){ + lprintf("MAIN",0,"Smear source to several Ns_sink\n"); + } if (mes_var.def_baryon) { lprintf("MAIN", 0, "Baryon masses\n"); @@ -330,6 +358,9 @@ int main(int argc, char *argv[]) { measure_spectrum_pt(tau, nm, m, mes_var.n_mom, i, mes_var.precision,DONTSTORE, NULL); } + if (mes_var.smearing_ss){ + measure_smearing_ss(mes_var.source_t,mes_var.source_x,mes_var.source_y,mes_var.source_z,nm,m,mes_var.n_mom,mes_var.nhits_2pt,i,mes_var.precision,mes_var.smear_epsilon_source,mes_var.smear_N_source,mes_var.smear_epsilon_sink,mes_var.smear_N_sink, mes_var.APE_epsilon, mes_var.APE_N); + } if (mes_var.def_baryon) { #if NG == 3 From f7b3468a8a22b867c56d4022fd4a7d7b05a1d145 Mon Sep 17 00:00:00 2001 From: Fabian Zierler Date: Mon, 27 May 2024 09:03:23 +0100 Subject: [PATCH 02/17] switch over to new complex numbers --- LibHR/Observables/avr_plaquette.c | 87 ++++++++++++++++-------------- LibHR/Observables/measure_mesons.c | 6 +-- LibHR/Observables/polyakov.c | 13 +++-- LibHR/Utils/suN_utils.c | 71 ++++-------------------- Spectrum/Makefile | 2 +- 5 files changed, 67 insertions(+), 112 deletions(-) diff --git a/LibHR/Observables/avr_plaquette.c b/LibHR/Observables/avr_plaquette.c index 1d8994c57..823b38283 100644 --- a/LibHR/Observables/avr_plaquette.c +++ b/LibHR/Observables/avr_plaquette.c @@ -450,6 +450,8 @@ void cplaq_APE(complex *ret,int ix,int mu,int nu) { int iy,iz; suNg *v1,*v2,*v3,*v4,w1,w2,w3; + double tmpre = 0.; + double tmpim = 0.; iy=iup(ix,mu); iz=iup(ix,nu); @@ -463,17 +465,16 @@ void cplaq_APE(complex *ret,int ix,int mu,int nu) _suNg_times_suNg(w2,(*v4),(*v3)); _suNg_times_suNg_dagger(w3,w1,w2); - _suNg_trace_re(ret->re,w3); + _suNg_trace_re(tmpre, w3); #ifdef GAUGE_SON - ret->im=0; + _suNg_trace_im(tmpim, w3); #else - _suNg_trace_im(ret->im,w3); + *ret = tmpre + I * tmpim; #endif #ifdef PLAQ_WEIGHTS if(plaq_weight!=NULL) { - ret->re *= plaq_weight[ix*16+mu*4+nu]; - ret->im *= plaq_weight[ix*16+mu*4+nu]; + *ret *= plaq_weight[ix * 16 + mu * 4 + nu];; } #endif @@ -513,13 +514,24 @@ double avr_spacial_plaquette_APE() void full_plaquette_APE() { - complex pa[6]; - double r0re = 0, r0im = 0; - double r1re = 0, r1im = 0; - double r2re = 0, r2im = 0; - double r3re = 0, r3im = 0; - double r4re = 0, r4im = 0; - double r5re = 0, r5im = 0; + + static double complex pa[6]; + static double complex r0; + static double complex r1; + static double complex r2; + static double complex r3; + static double complex r4; + static double complex r5; + + _OMP_PRAGMA(single) + { + r0 = 0.; + r1 = 0.; + r2 = 0.; + r3 = 0.; + r4 = 0.; + r5 = 0.; + } _PIECE_FOR(&glattice,ixp) { @@ -536,55 +548,50 @@ void full_plaquette_APE() complex tmp; cplaq_APE(&tmp,ix,1,0); - r0re += tmp.re; - r0im += tmp.im; + r0 += tmp; cplaq_APE(&tmp,ix,2,0); - r1re += tmp.re; - r1im += tmp.im; + r1 += tmp; cplaq_APE(&tmp,ix,2,1); - r2re += tmp.re; - r2im += tmp.im; + r2 += tmp; cplaq_APE(&tmp,ix,3,0); - r3re += tmp.re; - r3im += tmp.im; + r3 += tmp; cplaq_APE(&tmp,ix,3,1); - r4re += tmp.re; - r4im += tmp.im; + r4 += tmp; cplaq_APE(&tmp,ix,3,2); - r5re += tmp.re; - r5im += tmp.im; + r5 += tmp; } } - pa[0].re=r0re; pa[0].im=r0im; - pa[1].re=r1re; pa[1].im=r1im; - pa[2].re=r2re; pa[2].im=r2im; - pa[3].re=r3re; pa[3].im=r3im; - pa[4].re=r4re; pa[4].im=r4im; - pa[5].re=r5re; pa[5].im=r5im; + _OMP_PRAGMA(single) + { + pa[0] = r0; + pa[1] = r1; + pa[2] = r2; + pa[3] = r3; + pa[4] = r4; + pa[5] = r5; + } global_sum((double*)pa,12); for(int k = 0; k < 6; k++) { #ifdef BC_T_OPEN - pa[k].re /= NG*GLB_VOLUME*(GLB_T-1)/GLB_T; - pa[k].im /= NG*GLB_VOLUME*(GLB_T-1)/GLB_T; + pa[k] /= NG*GLB_VOLUME*(GLB_T-1)/GLB_T; #else - pa[k].re /= NG*GLB_VOLUME; - pa[k].im /= NG*GLB_VOLUME; + pa[k] /= NG*GLB_VOLUME; #endif } - lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",1,0,pa[0].re,pa[0].im); - lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",2,0,pa[1].re,pa[1].im); - lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",2,1,pa[2].re,pa[2].im); - lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",3,0,pa[3].re,pa[3].im); - lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",3,1,pa[4].re,pa[4].im); - lprintf("PLAQ",0,"Plaq(%d,%d) = ( %f , %f )\n",3,2,pa[5].re,pa[5].im); + lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 1, 0, creal(pa[0]), cimag(pa[0])); + lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 2, 0, creal(pa[1]), cimag(pa[1])); + lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 2, 1, creal(pa[2]), cimag(pa[2])); + lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 0, creal(pa[3]), cimag(pa[3])); + lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 1, creal(pa[4]), cimag(pa[4])); + lprintf("PLAQ", 0, "Plaq(%d,%d) = ( %f , %f )\n", 3, 2, creal(pa[5]), cimag(pa[5])); } \ No newline at end of file diff --git a/LibHR/Observables/measure_mesons.c b/LibHR/Observables/measure_mesons.c index 8cf6ee93f..8e0a97c93 100644 --- a/LibHR/Observables/measure_mesons.c +++ b/LibHR/Observables/measure_mesons.c @@ -690,8 +690,7 @@ void smeared_propagator(spinor_field* psi, int nm , double epsilon){ for (int a=0;ac[alpha].c[b].re = sp_smeared.c[b].c[alpha].c[beta].c[a].re; - _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b].im = sp_smeared.c[b].c[alpha].c[beta].c[a].im; + _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[b].c[alpha].c[beta].c[a]; } } @@ -775,8 +774,7 @@ void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon){ for (int a=0;ac[alpha].c[b].re = sp_smeared.c[b].c[alpha].c[beta].c[a].re; - _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b].im = sp_smeared.c[b].c[alpha].c[beta].c[a].im; + _FIELD_AT(&smeared_psi[a*4*nm+beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[b].c[alpha].c[beta].c[a]; } } diff --git a/LibHR/Observables/polyakov.c b/LibHR/Observables/polyakov.c index ae34d4075..095ba2de0 100644 --- a/LibHR/Observables/polyakov.c +++ b/LibHR/Observables/polyakov.c @@ -468,7 +468,7 @@ void polyakov_APE() { suNg *lp; suNg *bp; suNg tmp; - complex poly; + double complex poly; double adjpoly; double dtmp; #ifdef WITH_MPI @@ -672,7 +672,7 @@ void polyakov_APE() { /* trace and average */ - poly.re=poly.im=0.0; + _complex_0(poly); adjpoly=0.; i3d=0; @@ -681,11 +681,11 @@ void polyakov_APE() { for(x[(mu+3)%4]=0; x[(mu+3)%4]c[i].re = S.c[i].re; - g_out->c[i].im = S.c[i].im; + g_out->c[i] = S.c[i]; } } void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11){ - complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array + double complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array double UMAG; // [ 0 1 ] int i1,i2,i3,i4; // [ 2 3 ] --> [0, 1, 2, 3] @@ -543,7 +494,7 @@ void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11){ _complex_mul_star(ztmp1, F11, F11); _complex_mul_star(ztmp2, F12, F12); - UMAG = sqrt(ztmp1.re+ztmp2.re); + UMAG = sqrt(creal(ztmp1)+creal(ztmp2)); UMAG = 1./UMAG; _complex_mulr(F11, UMAG, F11); @@ -558,10 +509,10 @@ void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11){ vmxsu2(i1col,i2col,B11,A11); } -void vmxsu2(int i1, int i2, suNgfull* A, complex B[4]){ +void vmxsu2(int i1, int i2, suNgfull* A, double complex B[4]){ - complex C[NG*2]; - complex ztmp1,ztmp2; + double complex C[NG*2]; + double complex ztmp1,ztmp2; for (int i=0; ic[i + i1*NG], B[0]); @@ -581,7 +532,7 @@ void vmxsu2(int i1, int i2, suNgfull* A, complex B[4]){ void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11){ - complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array + double complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array double UMAG; // [ 0 1 ] int i1,i2,i3,i4; // [ 2 3 ] --> [0, 1, 2, 3] @@ -597,7 +548,7 @@ void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11){ _complex_mul_star(ztmp1, F11, F11); _complex_mul_star(ztmp2, F12, F12); - UMAG = sqrt(ztmp1.re+ztmp2.re); + UMAG = sqrt(creal(ztmp1)+creal(ztmp2)); UMAG = 1./UMAG; _complex_mulr(F11, UMAG, F11); @@ -615,7 +566,7 @@ void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11){ void vmxsu2_tau(int n1, int n2, suNgfull* A, complex B[4]){ suNgfull C; - complex ztmp1, ztmp2; + double complex ztmp1, ztmp2; _suNgfull_mul(C, 1., *A); diff --git a/Spectrum/Makefile b/Spectrum/Makefile index f4b03e703..c324409b1 100644 --- a/Spectrum/Makefile +++ b/Spectrum/Makefile @@ -2,7 +2,7 @@ TOPDIR = .. MKDIR = $(TOPDIR)/Make #EXES = random_cnfg random_spinor mk_mesons_with_z2semwall mk_mesons_with_z2semwall_new measure_spectrum measure_formfactor #mk_sfcoupling #trunc_mesons mk_mesons -EXES=measure_spectrum mk_sfcorrelators +EXES=measure_spectrum From 206a9a798dc9e2cf5d4c50b56d0443a2e478579b Mon Sep 17 00:00:00 2001 From: Fabian Zierler Date: Mon, 20 May 2024 18:13:57 +0100 Subject: [PATCH 03/17] rename source-and-sink smearing measurement code --- Include/observables.h | 2 +- LibHR/Observables/meson_measurements.c | 2 +- Spectrum/measure_spectrum.c | 10 +++++----- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Include/observables.h b/Include/observables.h index b986607f1..47da8823e 100644 --- a/Include/observables.h +++ b/Include/observables.h @@ -315,6 +315,6 @@ void smeared_propagator(spinor_field* psi, int nm , double epsilon); void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon); /* Measurements*/ -void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N); +void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N); #endif diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index ad23e21d7..987307b08 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -766,7 +766,7 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m /* Updates of Wuppertal smearing source in 2020*/ -void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N){ +void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N){ char label[20]; spinor_field* source = alloc_spinor_field_f(4,&glattice); diff --git a/Spectrum/measure_spectrum.c b/Spectrum/measure_spectrum.c index a8d5f5d2f..34bfd3406 100755 --- a/Spectrum/measure_spectrum.c +++ b/Spectrum/measure_spectrum.c @@ -78,7 +78,7 @@ typedef struct _input_mesons int source_x; int source_y; int source_z; - int smearing_ss; + int smearing_source_sink; double smear_epsilon_source; int smear_N_source; double smear_epsilon_sink; @@ -132,7 +132,7 @@ typedef struct _input_mesons {"Smearing Source x", "mes:source_x = %d",INT_T, &(varname).source_x}, \ {"Smearing Source y", "mes:source_y = %d",INT_T, &(varname).source_y}, \ {"Smearing Source z", "mes:source_z = %d",INT_T, &(varname).source_z}, \ - {"Smearing a lot at once", "mes:smearing_ss = %d",INT_T, &(varname).smearing_ss}, \ + {"Smearing a lot at once", "mes:smearing_source_sink = %d",INT_T, &(varname).smearing_source_sink}, \ {"smearing source step size", "mes:smear_epsilon_source = %lf", DOUBLE_T, &(varname).smear_epsilon_source}, \ {"Number of smeared source steps", "mes:smear_N_source = %d",INT_T, &(varname).smear_N_source}, \ {"smearing sink step size", "mes:smear_epsilon_sink = %lf", DOUBLE_T, &(varname).smear_epsilon_sink}, \ @@ -233,7 +233,7 @@ int main(int argc, char *argv[]) if (mes_var.APE_N != 0){ lprintf("MAIN",0,"APE smearing activate on Gaussian smearing\n"); } - if (mes_var.smearing_ss){ + if (mes_var.smearing_source_sink){ lprintf("MAIN",0,"Smear source to several Ns_sink\n"); } if (mes_var.def_baryon) @@ -358,8 +358,8 @@ int main(int argc, char *argv[]) { measure_spectrum_pt(tau, nm, m, mes_var.n_mom, i, mes_var.precision,DONTSTORE, NULL); } - if (mes_var.smearing_ss){ - measure_smearing_ss(mes_var.source_t,mes_var.source_x,mes_var.source_y,mes_var.source_z,nm,m,mes_var.n_mom,mes_var.nhits_2pt,i,mes_var.precision,mes_var.smear_epsilon_source,mes_var.smear_N_source,mes_var.smear_epsilon_sink,mes_var.smear_N_sink, mes_var.APE_epsilon, mes_var.APE_N); + if (mes_var.smearing_source_sink){ + measure_smearing_source_sink(mes_var.source_t,mes_var.source_x,mes_var.source_y,mes_var.source_z,nm,m,mes_var.n_mom,mes_var.nhits_2pt,i,mes_var.precision,mes_var.smear_epsilon_source,mes_var.smear_N_source,mes_var.smear_epsilon_sink,mes_var.smear_N_sink, mes_var.APE_epsilon, mes_var.APE_N); } if (mes_var.def_baryon) { From 133899ef26589349bd0aa7f95020930345a776d1 Mon Sep 17 00:00:00 2001 From: Fabian Zierler Date: Mon, 20 May 2024 18:30:20 +0100 Subject: [PATCH 04/17] perform all sink smearing levels at once --- Include/observables.h | 2 +- LibHR/Observables/meson_measurements.c | 16 +++++++++------- Spectrum/measure_spectrum.c | 6 ++++-- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/Include/observables.h b/Include/observables.h index 47da8823e..ff2433f2d 100644 --- a/Include/observables.h +++ b/Include/observables.h @@ -315,6 +315,6 @@ void smeared_propagator(spinor_field* psi, int nm , double epsilon); void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon); /* Measurements*/ -void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N); +void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff); #endif diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index 987307b08..4cf452fd0 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -766,19 +766,21 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m /* Updates of Wuppertal smearing source in 2020*/ -void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N){ +void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){ char label[20]; - spinor_field* source = alloc_spinor_field_f(4,&glattice); - spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice); - - init_propagator_eo(nm,m,precision); - if (APE_N != 0){ APE_smearing(APE_epsilon, APE_N); represent_gauge_field_APE(); } + int Nsmear_source; + for (Nsmear_source=0; Nsmear_source <= Nsmear_source_max; Nsmear_source += N_diff){ + + spinor_field* source = alloc_spinor_field_f(4,&glattice); + spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice); + init_propagator_eo(nm,m,precision); + int k; for (k=0;k Date: Mon, 20 May 2024 18:53:49 +0100 Subject: [PATCH 05/17] add smeared volume sources, smearing for disconnected pieces --- Include/observables.h | 7 +- Include/spectrum.h | 2 + LibHR/Observables/measure_mesons.c | 211 +++++++++++++++++++++++++ LibHR/Observables/meson_measurements.c | 202 +++++++++++++++++------ LibHR/Observables/sources.c | 174 ++++++++++++++++++++ Spectrum/measure_spectrum.c | 2 + 6 files changed, 551 insertions(+), 47 deletions(-) diff --git a/Include/observables.h b/Include/observables.h index ff2433f2d..14b417afe 100644 --- a/Include/observables.h +++ b/Include/observables.h @@ -306,13 +306,18 @@ void polyakov_APE(); /* source smearing*/ void smearing_function(spinor_field *source, int tau, int color, double epsilon); void smearing_function_with_APE(spinor_field *source, int tau, int color, double epsilon); +void smearing_function_volume(spinor_field *source, double epsilon); +void smearing_function_volume_with_APE(spinor_field *source, double epsilon); void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear); void create_smeared_source_with_APE(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear); +void create_noise_source_equal_eo_smeared_with_APE(spinor_field *source, double epsilon, int Nsmear); +void create_noise_source_equal_eo_smeared(spinor_field *source, double epsilon, int Nsmear); /* Sink smearing*/ - void smeared_propagator(spinor_field* psi, int nm , double epsilon); void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon); +void smeared_propagator_volume(spinor_field* psi, int nm , double epsilon); +void smeared_propagator_volume_with_APE(spinor_field* psi, int nm , double epsilon); /* Measurements*/ void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff); diff --git a/Include/spectrum.h b/Include/spectrum.h index 7a77cb2e0..d1021515a 100644 --- a/Include/spectrum.h +++ b/Include/spectrum.h @@ -19,6 +19,8 @@ void measure_spectrum_semwall(int nm, double* m, int nhits,int conf_num, double precision,storage_switch swc, data_storage_array **ret); void measure_spectrum_discon_semwall(int nm, double* m, int nhits,int conf_num, double precision,storage_switch swc, data_storage_array **ret); +void measure_spectrum_discon_semwall_smeared(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff); +void measure_spectrum_discon_semwall_smeared_single_inversion(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff); void measure_spectrum_discon_gfwall(int nm, double* m, int conf_num, double precision,storage_switch swc, data_storage_array **ret); void measure_spectrum_discon_volume(int nm, double* m, int conf_num, double precision, int dil,storage_switch swc, data_storage_array **ret); void measure_spectrum_gfwall(int nm, double* m, int conf_num, double precision,storage_switch swc, data_storage_array **ret); diff --git a/LibHR/Observables/measure_mesons.c b/LibHR/Observables/measure_mesons.c index 8e0a97c93..46ad08f50 100644 --- a/LibHR/Observables/measure_mesons.c +++ b/LibHR/Observables/measure_mesons.c @@ -790,4 +790,215 @@ void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon){ } free_spinor_field_f(smeared_psi); } +// FZ (2024): Smeared propagator with a single source +void smeared_propagator_volume(spinor_field* psi, int nm , double epsilon){ + int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; + double norm_factor = 1./(1.+6.*epsilon); + suNf_spin_matrix sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down; + _spinmatrix_zero(sp_smeared); + + spinor_field* smeared_psi = alloc_spinor_field_f(4*nm,&glattice); + + for(i=0; ic[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b]; + _FIELD_AT(&smeared_psi[beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b]; + } + } + + } + } + } + for (int i=0;i<4*nm;i++){ + spinor_field_copy_f(psi + i, smeared_psi + i); + } + for (int i=0;i<4*nm;i++){ + start_sf_sendrecv(psi + i); + complete_sf_sendrecv(psi + i); + } + free_spinor_field_f(smeared_psi); +} +void smeared_propagator_volume_with_APE(spinor_field* psi, int nm , double epsilon){ + int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; + double norm_factor = 1./(1.+6.*epsilon); + suNf_spin_matrix sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down; + _spinmatrix_zero(sp_smeared); + + spinor_field* smeared_psi = alloc_spinor_field_f(4*nm,&glattice); + + for(i=0; ic[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b]; + _FIELD_AT(&smeared_psi[beta*nm+i],ix)->c[alpha].c[b] = sp_smeared.c[beta].c[alpha].c[b]; + } + } + + } + } + } + for (int i=0;i<4*nm;i++){ + spinor_field_copy_f(psi + i, smeared_psi + i); + } + for (int i=0;i<4*nm;i++){ + start_sf_sendrecv(psi + i); + complete_sf_sendrecv(psi + i); + } + free_spinor_field_f(smeared_psi); +} diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index 4cf452fd0..5a5e81f16 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -768,7 +768,7 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){ - char label[20]; + char label[20]; if (APE_N != 0){ APE_smearing(APE_epsilon, APE_N); represent_gauge_field_APE(); @@ -781,58 +781,168 @@ void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice); init_propagator_eo(nm,m,precision); - int k; - for (k=0;k1){ - measure_point_mesons_momenta(meson_correlators,prop+4*k, source, nm, t, n_mom); - } - else{ - measure_mesons(meson_correlators,prop+4*k, source, nm, t); - } + int k; + for (k=0;k1){ + measure_point_mesons_momenta(meson_correlators,prop+4*k, source, nm, t, n_mom); + } + else{ + measure_mesons(meson_correlators,prop+4*k, source, nm, t); + } + } + + sprintf(label, "source_N%d_sink_N%d",Nsmear_source, 0); + print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); + + for (int N = 0; N1){ + measure_point_mesons_momenta(meson_correlators,prop+4*c, source, nm, t, n_mom); + } + else{ + measure_mesons(meson_correlators,prop+4*c, source, nm, t); + } + } + sprintf(label, "source_N%d_sink_N%d",Nsmear_source, N+N_diff); + print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); + } + + free_propagator_eo(); + free_spinor_field_f(source); + free_spinor_field_f(prop); } +} + + +// So far, only implement the measurement for mesons at rest +// The sources are not localised on any lattice site but spread out +/* Smearing for disconnected pieces for singlets (FZ Jan. 2024)*/ +void measure_spectrum_discon_semwall_smeared_single_inversion(int nm, double* m, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double APE_epsilon, int APE_N, int N_diff){ - sprintf(label, "source_N%d_sink_N%d",Nsmear_source, 0); - print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); + char label[32]; + int k,beta; + spinor_field* source = alloc_spinor_field_f(4,&glattice); + spinor_field* prop = alloc_spinor_field_f(4*nm, &glattice); - for (int N = 0; N1){ - measure_point_mesons_momenta(meson_correlators,prop+4*c, source, nm, t, n_mom); - } - else{ - measure_mesons(meson_correlators,prop+4*c, source, nm, t); - } + init_propagator_eo(nm,m,precision); + if (APE_N != 0){ + APE_smearing(APE_epsilon, APE_N); + represent_gauge_field_APE(); + } + + lprintf("SMEAR",0,"source smearing with steps %d up to %d\n ",N_diff,Nsmear_source); + + for (k=0;kc[spin] = spinor_smeared.c[spin]; + } + } + + for (int beta=0;beta<4;++beta){ + spinor_field_copy_f(source + beta, smeared_source +beta); + } + + + for (int beta=0;beta<4;++beta){ + start_sf_sendrecv(source + beta); + complete_sf_sendrecv(source + beta); + } + free_spinor_field_f(smeared_source); +} + +void smearing_function_volume_with_APE(spinor_field *source, double epsilon){ + + int ix, t, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; + double norm_factor = 1./(1.+6.*epsilon); + suNf_spinor spinor_OG, spinor_smeared, spinor_tmp, spinor_right, spinor_left, spinor_front, spinor_back, spinor_up, spinor_down; + _spinor_zero_f(spinor_OG);_spinor_zero_f(spinor_tmp);_spinor_zero_f(spinor_right);_spinor_zero_f(spinor_left); + _spinor_zero_f(spinor_front);_spinor_zero_f(spinor_back);_spinor_zero_f(spinor_up);_spinor_zero_f(spinor_down); + + spinor_field* smeared_source = alloc_spinor_field_f(4, &glattice); + for (int beta=0;beta<4;++beta){ + spinor_field_zero_f(&smeared_source[beta]); + } + + represent_gauge_field_APE(); + + for (t=0; t Date: Mon, 20 May 2024 18:58:30 +0100 Subject: [PATCH 06/17] add example input file for meson measurements with smearing --- Spectrum/input_file_smeared | 74 +++++++++++++++++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 Spectrum/input_file_smeared diff --git a/Spectrum/input_file_smeared b/Spectrum/input_file_smeared new file mode 100644 index 000000000..3e8776e62 --- /dev/null +++ b/Spectrum/input_file_smeared @@ -0,0 +1,74 @@ +// Global variables +GLB_T = 8 +GLB_X = 4 +GLB_Y = 4 +GLB_Z = 4 +NP_T = 2 +NP_X = 1 +NP_Y = 1 +NP_Z = 1 + +// Replicas +N_REP = 1 + +// Random generator +rlx_level = 1 +rlx_seed = 11228 +rlx_start = new +rlx_state = ./out/rlx_state + +log:default = -1 +log:inverter = -1 +log:forcestat = 0 + +mes:precision = 1.e-14 +mes:nhits = 1 +mes:nhits_disc = 4 +mes:masses = -0.85 +mes:use_input_mass = 0 +mes:meas_mixed = 0 +mes:nhits_2pt = 1 +mes:momentum = 0 +mes:def_semwall = 0 +mes:def_semwall_nondeg = 0 +mes:def_point = 0 +mes:def_baryon = 0 +mes:def_gfwall = 0 +mes:ext_semwall = 0 +mes:ext_point = 0 +mes:dirichlet_semwall = 0 +mes:dirichlet_point = 0 +mes:dirichlet_gfwall = 0 +mes:discon_semwall = 0 +mes:discon_gfwall = 0 +mes:discon_volume = 0 +mes:dilution = 1 +mes:dirichlet_dt = 2 +mes:configlist = listfile_spectrum + +// Measure smeared meson spectrum +// enable source and sink smearing +mes:smearing_source_sink = 1 +mes:discon_semwall_smeared = 1 + +// point source location for source ans sink smearing +mes:source_t = 0 +mes:source_x = 0 +mes:source_y = 0 +mes:source_z = 0 + +// source and sink smearing parameters +mes:smear_epsilon_source = 0.2 +mes:smear_N_source = 50 +mes:smear_epsilon_sink = 0.2 +mes:smear_N_sink = 50 +mes:smear_N_step = 10 +//APE smearing activate if APE_N != 0 +mes:APE_epsilon = 0.4 +mes:APE_N = 10 + +//Whether to include the four fermion interaction +ff:on = false +//Parameters for hopping parameter expansion, 0 in degree_hopping to disable +mes:degree_hopping = 0 +mes:nhits_hopping = 5 From 1b204b666a23c4f18047cc233941dbdf8c8a3fcd Mon Sep 17 00:00:00 2001 From: Fabian Zierler Date: Sat, 1 Jun 2024 12:26:00 +0100 Subject: [PATCH 07/17] add missing call to smeared disconnected measurement --- Spectrum/measure_spectrum.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Spectrum/measure_spectrum.c b/Spectrum/measure_spectrum.c index f6cc50de2..e7a48b71a 100755 --- a/Spectrum/measure_spectrum.c +++ b/Spectrum/measure_spectrum.c @@ -405,6 +405,9 @@ int main(int argc, char *argv[]) { measure_spectrum_discon_semwall(nm, m, mes_var.nhits_disc, i, mes_var.precision,DONTSTORE, NULL); } + if (mes_var.discon_semwall_smeared){ + measure_spectrum_discon_semwall_smeared(nm,m,mes_var.nhits_disc,i,mes_var.precision,mes_var.smear_epsilon_source,mes_var.smear_N_source, mes_var.APE_epsilon, mes_var.APE_N, mes_var.smear_N_step); + } if (mes_var.discon_gfwall) { // measure_spectrum_discon_gfwall(nm,m,i,mes_var.precision,DONTSTORE, NULL); From 160e19e89d621cbb7317c25920fb48c18bbe9c58 Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:18:34 +0900 Subject: [PATCH 08/17] add comments for Wuppertal smearing and helper function 1. write descriptions for the function 2. define helper functions to increase readability 3. remove unused arguments `int color` in `smearing_function` --- Include/observables.h | 8 +- LibHR/Observables/meson_measurements.c | 134 +++++++++++++++---------- LibHR/Observables/sources.c | 59 +++++++++-- 3 files changed, 135 insertions(+), 66 deletions(-) diff --git a/Include/observables.h b/Include/observables.h index 14b417afe..624731199 100644 --- a/Include/observables.h +++ b/Include/observables.h @@ -304,8 +304,8 @@ void full_plaquette_APE(); void polyakov_APE(); /* source smearing*/ -void smearing_function(spinor_field *source, int tau, int color, double epsilon); -void smearing_function_with_APE(spinor_field *source, int tau, int color, double epsilon); +void smearing_function(spinor_field *source, int tau, double epsilon); +void smearing_function_with_APE(spinor_field *source, int tau, double epsilon); void smearing_function_volume(spinor_field *source, double epsilon); void smearing_function_volume_with_APE(spinor_field *source, double epsilon); void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear); @@ -320,6 +320,8 @@ void smeared_propagator_volume(spinor_field* psi, int nm , double epsilon); void smeared_propagator_volume_with_APE(spinor_field* psi, int nm , double epsilon); /* Measurements*/ -void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff); +void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff); +void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff); + #endif diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index 5a5e81f16..c1e056716 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -764,7 +764,62 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m -/* Updates of Wuppertal smearing source in 2020*/ +/* Updates for Wuppertal smearing source by HH in 2020*/ + +void create_smear_source_and_measure(spinor_field* source, spinor_field* prop, int t, int x, int y, int z, int nm, double* m, int n_mom, double epsilon_source, int Nsmear_source, int APE_N,int conf_num, char* label) { + /** + * This function creates a smeared source and measures the meson correlators. + * The APE smearing on the gauge links are applied if `APE_N` is not 0. + */ + for (int k = 0;k1){ + measure_point_mesons_momenta(meson_correlators,prop+4*k, source, nm, t, n_mom); + } else{ + measure_mesons(meson_correlators,prop+4*k, source, nm, t); + } + } + + sprintf(label, "source_N%d_sink_N%d",Nsmear_source, 0); + print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); + +} + + +void smear_sink_and_measure(spinor_field* prop, spinor_field* source, int nm, double* m, int t, int n_mom, int Nsmear_source, double epsilon_sink, int Nsmear_sink, int N_diff, int APE_N, int conf_num, char* label) { + /** + * This function smears the sink and measures the meson correlators. + * The APE smearing on the gauge links are applied if `APE_N` is not 0. + * The numbers of the smearing iterations at the sink are from 0 to `Nsmear_sink` with a step of `N_diff`. + */ + for (int N = 0; N1){ + measure_point_mesons_momenta(meson_correlators,prop+4*c, source, nm, t, n_mom); + } else { + measure_mesons(meson_correlators,prop+4*c, source, nm, t); + } + } + sprintf(label, "source_N%d_sink_N%d",Nsmear_source, N+N_diff); + print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); + } +} + void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){ @@ -774,62 +829,15 @@ void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, represent_gauge_field_APE(); } - int Nsmear_source; - for (Nsmear_source=0; Nsmear_source <= Nsmear_source_max; Nsmear_source += N_diff){ + for (int Nsmear_source=0; Nsmear_source <= Nsmear_source_max; Nsmear_source += N_diff){ spinor_field* source = alloc_spinor_field_f(4,&glattice); spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice); init_propagator_eo(nm,m,precision); - int k; - for (k=0;k1){ - measure_point_mesons_momenta(meson_correlators,prop+4*k, source, nm, t, n_mom); - } - else{ - measure_mesons(meson_correlators,prop+4*k, source, nm, t); - } - } - - sprintf(label, "source_N%d_sink_N%d",Nsmear_source, 0); - print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); - - for (int N = 0; N1){ - measure_point_mesons_momenta(meson_correlators,prop+4*c, source, nm, t, n_mom); - } - else{ - measure_mesons(meson_correlators,prop+4*c, source, nm, t); - } - } - sprintf(label, "source_N%d_sink_N%d",Nsmear_source, N+N_diff); - print_mesons(meson_correlators,1.,conf_num,nm,m,GLB_T,n_mom,label); - } + create_smear_source_and_measure(source, prop, t, x, y, z, nm, m, n_mom, epsilon_source, Nsmear_source, APE_N, conf_num, label); + + smear_sink_and_measure(prop, source, nm, m, t, n_mom, Nsmear_source, epsilon_sink, Nsmear_sink, N_diff, APE_N, conf_num, label); free_propagator_eo(); free_spinor_field_f(source); @@ -837,6 +845,28 @@ void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, } } +void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){ + + char label[20]; + if (APE_N != 0){ + APE_smearing(APE_epsilon, APE_N); + represent_gauge_field_APE(); + } + + spinor_field* source = alloc_spinor_field_f(4,&glattice); + spinor_field* prop = alloc_spinor_field_f(4*nm*NF, &glattice); + init_propagator_eo(nm,m,precision); + + create_smear_source_and_measure(source, prop, t, x, y, z, nm, m, n_mom, epsilon_source, Nsmear_source, APE_N, conf_num, label); + + smear_sink_and_measure(prop, source, nm, m, t, n_mom, Nsmear_source, epsilon_sink, Nsmear_sink, N_diff, APE_N, conf_num, label); + + free_propagator_eo(); + free_spinor_field_f(source); + free_spinor_field_f(prop); + +} + // So far, only implement the measurement for mesons at rest // The sources are not localised on any lattice site but spread out diff --git a/LibHR/Observables/sources.c b/LibHR/Observables/sources.c index 4dadbd65b..67fc0ddd3 100755 --- a/LibHR/Observables/sources.c +++ b/LibHR/Observables/sources.c @@ -89,6 +89,11 @@ static int random_tau() source[spin](t,x) = \sum_{x%p == 0} \delta_{s spin} 1_color z2_volume_source: source[spin](t,x) = Z(2) x Z(2) (no dilution) + smeared_source: + Wuppertal smeared source, see Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. + The involving gauge field is the original one. + smeared_source_with_APE: + Wuppertal smearing with APE smeared gauge field. \***************************************************************************/ void create_point_source(spinor_field *source, int tau, int color) { @@ -714,9 +719,20 @@ void zero_even_or_odd_site_spinorfield(spinor_field *source, int nspinor, int eo } -/* Updates of Wuppertal smearing source in 2020 */ - -void smearing_function(spinor_field *source, int tau, int color, double epsilon){ +/* Updates for Wuppertal smearing source by HH in 2020*/ + +void smearing_function(spinor_field *source, int tau, double epsilon){ + /** + * This function applies the Wuppertal smearing to the source field. + * The smearing is done in the spatial directions only, basing on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf: + * \Phi q(x) = 1/(1+2d\epsilon)[ q(x) + \epsilon\sum_{\mu=1}^{3} U_{\mu}(x)q(x+\mu) ] + * \Phi is the smearing function, q(x) is the source field, U_{\mu}(x) is the link variable in the \mu direction. + * d = 3 for spatial directions. + * + * q -> `source` + * \epsilon -> `epsilon` + * `tau` is the time slice where the source is located. + */ int ix, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; double norm_factor = 1./(1.+6.*epsilon); @@ -729,9 +745,10 @@ void smearing_function(spinor_field *source, int tau, int color, double epsilon) spinor_field_zero_f(&smeared_source[beta]); } - if(COORD[0]==tau/T){ + if(COORD[0]==tau/T){ // apply the smearing only on the time slice where the source is located. for (x=0; x `pu_gauge_APE_f` + */ int ix, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; double norm_factor = 1./(1.+6.*epsilon); @@ -861,6 +885,13 @@ void smearing_function_with_APE(spinor_field *source, int tau, int color, double } void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear){ + /** + * This function creates a smeared source field at a given position (t,x,y,z) with a given color by applying + * the Wuppertal smearing function based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. + * `Nsmear` is the iteration number at the source + * `epsilon` is the step size. + * The gauge field is the original one. + */ int beta; @@ -875,12 +906,18 @@ void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int for (int n=0;n Date: Mon, 6 Jan 2025 15:18:43 +0900 Subject: [PATCH 09/17] Update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 6257730b3..bbadc0563 100644 --- a/.gitignore +++ b/.gitignore @@ -28,3 +28,4 @@ Hirep-git.code-workspace **/clangd .ycm_extra_conf.py* compile_commands.json +Spectrum/measure_spectrum From 2e987b3b6a6b7d2d9d117b09dfc22fc7ff012851 Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:53:42 +0900 Subject: [PATCH 10/17] refine comments --- LibHR/Observables/sources.c | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/LibHR/Observables/sources.c b/LibHR/Observables/sources.c index 67fc0ddd3..85cc31daa 100755 --- a/LibHR/Observables/sources.c +++ b/LibHR/Observables/sources.c @@ -729,9 +729,9 @@ void smearing_function(spinor_field *source, int tau, double epsilon){ * \Phi is the smearing function, q(x) is the source field, U_{\mu}(x) is the link variable in the \mu direction. * d = 3 for spatial directions. * - * q -> `source` - * \epsilon -> `epsilon` - * `tau` is the time slice where the source is located. + * @param source -> the source field, q(x). + * @param epsilon -> the smearing parameter, \epsilon + * @param tau is the time slice where the source is located. */ int ix, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; @@ -888,8 +888,9 @@ void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int /** * This function creates a smeared source field at a given position (t,x,y,z) with a given color by applying * the Wuppertal smearing function based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. - * `Nsmear` is the iteration number at the source - * `epsilon` is the step size. + * @param Nsmear is the iteration number at the source + * @param epsilon is the step size. + * * The gauge field is the original one. */ @@ -914,8 +915,10 @@ void create_smeared_source_with_APE(spinor_field *source, int t, int x, int y, i /** * This function creates a smeared source field at a given position (t,x,y,z) with a given color by applying * the Wuppertal smearing function based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. - * `Nsmear` is the iteration number at the source - * `epsilon` is the step size. + * + * @param Nsmear is the iteration number at the source + * @param epsilon is the step size. + * * The gauge field is the APE smeared one. */ int beta; From 0588c22111ebbce4616d7762fc9fc0e049b9f919 Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Mon, 6 Jan 2025 15:53:52 +0900 Subject: [PATCH 11/17] add comments --- LibHR/Utils/APE_smearing.c | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/LibHR/Utils/APE_smearing.c b/LibHR/Utils/APE_smearing.c index 91203c5e1..cec4d5bfe 100644 --- a/LibHR/Utils/APE_smearing.c +++ b/LibHR/Utils/APE_smearing.c @@ -1,9 +1,9 @@ // // APE_smearing.c -// hr // -// Created by Paul Xiao on 2020/11/2. -// Copyright © 2020 Claudio Pica. All rights reserved. +// Code added by HH on 2020.11.2 for APE smearing, +// based on the function written in Fortran provided by Davide Vadacchino. +// // #include @@ -23,14 +23,21 @@ #include "update.h" void APE_smearing(double smear_val, int Nsmear){ + /** + * APE smearing function is used to smooth the gauge field, we take the defination + * from Eq. (16) in https://arxiv.org/abs/hep-lat/0409141v2. + * + * @param smear_val: Smearing parameter \alpha in the reference + * @param Nsmear: Number of smearing iterations. + */ suNg staple, tr1, tr2; int ix, iy, iz, it; int mu, nu, mid, midpmu, midpnu, midmnu, midpmumnu; suNg v, vout, vtmp; - double sm1 = 1. - smear_val; - double sm2 = smear_val / 6. ; + double sm1 = 1. - smear_val; // ( 1 - \alpha ) + double sm2 = smear_val / 6. ; // ( \alpha \ 6 ) suNg_field *u_gauge_tmp = alloc_gfield(&glattice); u_gauge_APE = alloc_gfield(&glattice); @@ -91,7 +98,7 @@ void APE_smearing(double smear_val, int Nsmear){ _suNg_add_assign(vout,tr2); _suNg_mul(vtmp, 1., vout); - project_to_suNg(&vout); + project_to_suNg(&vout); // Project to the group manifold #ifdef GAUGE_SPN cooling_SPN(&v, &vout, &vtmp, 5); // max Re Tr(U V^+) From f5f55a6d89291a58fb5713927f13b831fb8f367a Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Mon, 6 Jan 2025 16:11:09 +0900 Subject: [PATCH 12/17] refine comments --- LibHR/Observables/measure_mesons.c | 10 ++++++++++ LibHR/Observables/meson_measurements.c | 20 ++++++++++++++++---- LibHR/Observables/sources.c | 12 ++++++------ LibHR/Utils/APE_smearing.c | 2 +- 4 files changed, 33 insertions(+), 11 deletions(-) diff --git a/LibHR/Observables/measure_mesons.c b/LibHR/Observables/measure_mesons.c index 46ad08f50..c7d5bdf2a 100644 --- a/LibHR/Observables/measure_mesons.c +++ b/LibHR/Observables/measure_mesons.c @@ -628,6 +628,16 @@ void print_mesons(meson_observable *mo, double norm, int conf, int nm, double *m /* Wuppertal smearing*/ void smeared_propagator(spinor_field* psi, int nm , double epsilon){ + /** + * @brief Smears the propagator field `psi` with the Wuppertal smearing with parameter `epsilon` and one iteration. + * The smearing function is defined as Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. + * See also disucssion in the paragraph below Eq. (21) in the same paper. + * + * @param psi The propagator field to be smeared, $S$ in the paper. + * @param nm The number of masses + * @param epsilon The smearing step size, $\epsilon$ in Eq. (9). + */ + int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; double norm_factor = 1./(1.+6.*epsilon); suNf_propagator sp_OG, sp_smeared, sp_tmp, sp_right, sp_left, sp_front, sp_back, sp_up, sp_down; diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index c1e056716..e42bf4001 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -764,7 +764,7 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m -/* Updates for Wuppertal smearing source by HH in 2020*/ +/* Updates for Wuppertal smearing by HH in 2020*/ void create_smear_source_and_measure(spinor_field* source, spinor_field* prop, int t, int x, int y, int z, int nm, double* m, int n_mom, double epsilon_source, int Nsmear_source, int APE_N,int conf_num, char* label) { /** @@ -793,7 +793,7 @@ void create_smear_source_and_measure(spinor_field* source, spinor_field* prop, i void smear_sink_and_measure(spinor_field* prop, spinor_field* source, int nm, double* m, int t, int n_mom, int Nsmear_source, double epsilon_sink, int Nsmear_sink, int N_diff, int APE_N, int conf_num, char* label) { /** - * This function smears the sink and measures the meson correlators. + * This function smears the sink and measures the meson correlators, with the smearing steps from 0 to `Nsmear_sink` with a step of `N_diff`. * The APE smearing on the gauge links are applied if `APE_N` is not 0. * The numbers of the smearing iterations at the sink are from 0 to `Nsmear_sink` with a step of `N_diff`. */ @@ -822,7 +822,12 @@ void smear_sink_and_measure(spinor_field* prop, spinor_field* source, int nm, do void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source_max, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){ - + /** + * This function measures the meson correlators with smeared source and sink. + * It measures the meson correlators for all combinations of the source and sink smearing steps, [0, Nsmear_source_max] and [0, Nsmear_sink], respectively. + * + * The APE smearing on the gauge links are applied if `APE_N` is not 0. + */ char label[20]; if (APE_N != 0){ APE_smearing(APE_epsilon, APE_N); @@ -846,7 +851,14 @@ void measure_smearing_source_sink(int t, int x, int y, int z, int nm, double* m, } void measure_smearing_ss(int t, int x, int y, int z, int nm, double* m, int n_mom, int nhits, int conf_num, double precision, double epsilon_source, int Nsmear_source, double epsilon_sink, int Nsmear_sink, double APE_epsilon, int APE_N, int N_diff){ - + /** + * This function measures the meson correlators with smeared source and sink. + * It measures the meson correlators for a specifc source smearing steps, Nsmear_source, + * while the sink smearing step is [0, Nsmear_sink]. + * + * The APE smearing on the gauge links are applied if `APE_N` is not 0. + */ + char label[20]; if (APE_N != 0){ APE_smearing(APE_epsilon, APE_N); diff --git a/LibHR/Observables/sources.c b/LibHR/Observables/sources.c index 85cc31daa..68c6ee8d7 100755 --- a/LibHR/Observables/sources.c +++ b/LibHR/Observables/sources.c @@ -719,18 +719,18 @@ void zero_even_or_odd_site_spinorfield(spinor_field *source, int nspinor, int eo } -/* Updates for Wuppertal smearing source by HH in 2020*/ +/* Updates for Wuppertal smearing by HH in 2020*/ void smearing_function(spinor_field *source, int tau, double epsilon){ /** * This function applies the Wuppertal smearing to the source field. * The smearing is done in the spatial directions only, basing on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf: - * \Phi q(x) = 1/(1+2d\epsilon)[ q(x) + \epsilon\sum_{\mu=1}^{3} U_{\mu}(x)q(x+\mu) ] - * \Phi is the smearing function, q(x) is the source field, U_{\mu}(x) is the link variable in the \mu direction. - * d = 3 for spatial directions. + * $\Phi q(x) = 1/(1+2d\epsilon)[ q(x) + \epsilon\sum_{\mu=1}^{3} U_{\mu}(x)q(x+\mu) ]$ + * $\Phi$ is the smearing function, $q(x)$ is the source field, $U_{\mu}(x)$ is the link variable in the $\mu$ direction. + * $d = 3$ for spatial directions. * - * @param source -> the source field, q(x). - * @param epsilon -> the smearing parameter, \epsilon + * @param source -> the source field, $q(x)$. + * @param epsilon -> the smearing parameter, $\epsilon$. * @param tau is the time slice where the source is located. */ diff --git a/LibHR/Utils/APE_smearing.c b/LibHR/Utils/APE_smearing.c index cec4d5bfe..2d63cdfcb 100644 --- a/LibHR/Utils/APE_smearing.c +++ b/LibHR/Utils/APE_smearing.c @@ -2,7 +2,7 @@ // APE_smearing.c // // Code added by HH on 2020.11.2 for APE smearing, -// based on the function written in Fortran provided by Davide Vadacchino. +// based on the functions written in Fortran provided by Davide Vadacchino. // // From ab43caafd5abc2950a4839b15677de9d3af672e0 Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Mon, 6 Jan 2025 16:26:03 +0900 Subject: [PATCH 13/17] refine comments --- LibHR/Observables/measure_mesons.c | 16 +++++++++++++++- LibHR/Observables/meson_measurements.c | 7 ++++--- LibHR/Utils/APE_smearing.c | 2 +- 3 files changed, 20 insertions(+), 5 deletions(-) diff --git a/LibHR/Observables/measure_mesons.c b/LibHR/Observables/measure_mesons.c index c7d5bdf2a..d3e2c0108 100644 --- a/LibHR/Observables/measure_mesons.c +++ b/LibHR/Observables/measure_mesons.c @@ -629,7 +629,7 @@ void print_mesons(meson_observable *mo, double norm, int conf, int nm, double *m void smeared_propagator(spinor_field* psi, int nm , double epsilon){ /** - * @brief Smears the propagator field `psi` with the Wuppertal smearing with parameter `epsilon` and one iteration. + * @brief Smears the propagator field `psi` with the Wuppertal smearingwith parameter `epsilon` and one iteration, without APE smearing. * The smearing function is defined as Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. * See also disucssion in the paragraph below Eq. (21) in the same paper. * @@ -719,6 +719,18 @@ void smeared_propagator(spinor_field* psi, int nm , double epsilon){ } void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon){ + /** + * @brief Smears the propagator field `psi` with the Wuppertal smearing with parameter `epsilon` and one iteration. + * The smearing function is defined as Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. + * See also disucssion in the paragraph below Eq. (21) in the same paper. + * + * @param psi The propagator field to be smeared, $S$ in the paper. + * @param nm The number of masses + * @param epsilon The smearing step size, $\epsilon$ in Eq. (9). + * + * The gauge field is smeared with APE smearing before the Wuppertal smearing of the propagator field. + */ + int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; double norm_factor = 1./(1.+6.*epsilon); @@ -800,6 +812,8 @@ void smeared_propagator_with_APE(spinor_field* psi, int nm , double epsilon){ } free_spinor_field_f(smeared_psi); } + + // FZ (2024): Smeared propagator with a single source void smeared_propagator_volume(spinor_field* psi, int nm , double epsilon){ int i,t, x, y, z, ix, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index e42bf4001..ea9c20263 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -768,7 +768,7 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m void create_smear_source_and_measure(spinor_field* source, spinor_field* prop, int t, int x, int y, int z, int nm, double* m, int n_mom, double epsilon_source, int Nsmear_source, int APE_N,int conf_num, char* label) { /** - * This function creates a smeared source and measures the meson correlators. + * This function creates a smeared source and measures the meson correlators with point sink. * The APE smearing on the gauge links are applied if `APE_N` is not 0. */ for (int k = 0;k Date: Thu, 9 Jan 2025 17:48:11 +0900 Subject: [PATCH 14/17] correct typo --- Include/observables.h | 2 +- LibHR/Utils/APE_smearing.c | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Include/observables.h b/Include/observables.h index 624731199..2395b4971 100644 --- a/Include/observables.h +++ b/Include/observables.h @@ -299,7 +299,7 @@ void ff_observables(); //APE smearing double plaq_APE(int ix,int mu,int nu); void cplaq_APE(complex *ret,int ix,int mu,int nu); -double avr_spacial_plaquette_APE(); +double avr_spatial_plaquette_APE(); void full_plaquette_APE(); void polyakov_APE(); diff --git a/LibHR/Utils/APE_smearing.c b/LibHR/Utils/APE_smearing.c index f70263105..35fbf932f 100644 --- a/LibHR/Utils/APE_smearing.c +++ b/LibHR/Utils/APE_smearing.c @@ -60,7 +60,7 @@ void APE_smearing(double smear_val, int Nsmear){ represent_gauge_field_APE(); lprintf("APE",0,"APE smearing with val=%f \n", smear_val); - lprintf("APE",0,"N=0

_s = %1.6f\n",avr_spacial_plaquette_APE()); + lprintf("APE",0,"N=0

_s = %1.6f\n",avr_spatial_plaquette_APE()); for (int N=0; N_s = %1.6f\n", Nsmear, avr_spacial_plaquette_APE()); + lprintf("APE",0,"N=%d

_s = %1.6f\n", Nsmear, avr_spatial_plaquette_APE()); polyakov_APE(); full_plaquette_APE(); lprintf("APE",0,"APE smearing END \n"); From de20a96d86bb7de4e4ffe3ddf5435ce487deb899 Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Thu, 9 Jan 2025 17:49:38 +0900 Subject: [PATCH 15/17] add comments --- LibHR/Observables/avr_plaquette.c | 48 +++++++++++++++++++++++++------ 1 file changed, 40 insertions(+), 8 deletions(-) diff --git a/LibHR/Observables/avr_plaquette.c b/LibHR/Observables/avr_plaquette.c index 823b38283..a6e9ba4f1 100644 --- a/LibHR/Observables/avr_plaquette.c +++ b/LibHR/Observables/avr_plaquette.c @@ -418,8 +418,17 @@ double complex avr_plaquette_wrk() } -double plaq_APE(int ix,int mu,int nu) -{ +double plaq_APE(int ix,int mu,int nu){ + /** + * @brief Computes the value of the plaquette at the given position in the mu-nu plane and return only the real part. + * Modified from `plaq(int ix, int mu, int nu)` by changing `pu_gauge` to `pu_gauge_APE` + * + * @param ix The position of the starting site. + * @param mu The first direction of the plaquette. + * @param nu The second direction of the plaquette. + * + * @return The real part of the computed plaquette value, possibly weighted if PLAQ_WEIGHTS is defined. + */ int iy,iz; double p; suNg *v1,*v2,*v3,*v4,w1,w2,w3; @@ -446,8 +455,19 @@ double plaq_APE(int ix,int mu,int nu) #endif } -void cplaq_APE(complex *ret,int ix,int mu,int nu) -{ + +void cplaq_APE(complex *ret,int ix,int mu,int nu){ + /** + * @brief Computes the value of the plaquette at the given position in the mu-nu plane. + * Modified from `cplaq(complex *ret,int ix,int mu,int nu)` by changing `pu_gauge` to `pu_gauge_APE` + * + * @param ret A pointer to store the resulting complex plaquette value. + * @param ix The position of the starting site. + * @param mu The first direction of the plaquette. + * @param nu The second direction of the plaquette. + * + * @note If PLAQ_WEIGHTS is defined, the result may include a weighting factor. + */ int iy,iz; suNg *v1,*v2,*v3,*v4,w1,w2,w3; double tmpre = 0.; @@ -481,8 +501,13 @@ void cplaq_APE(complex *ret,int ix,int mu,int nu) } -double avr_spacial_plaquette_APE() -{ +double avr_spatial_plaquette_APE(){ + /** + * @brief Computes the average (over the whole lattice) plaquette value over all spatial directions + * for the APE-smeared gauge links. + * Modified from `avr_spatial_plaquette()` by changing `plaq()` to `plaq_APE()` + * + */ double pa=0.; _PIECE_FOR(&glattice,ixp) { @@ -512,8 +537,15 @@ double avr_spacial_plaquette_APE() } -void full_plaquette_APE() -{ +void full_plaquette_APE(){ + /** + * @brief Computes and prints the average plaquette value over all directions, + * including the temporal direction, for the APE-smeared gauge links. + * Modified from `ull_plaquette())` by changing `cplaq(...)` to `cplaq_APE(...)` + * + * @note APE smearing is not applied in the temporal direction. + * As a result, plaquette values involving the temporal direction are zero. + */ static double complex pa[6]; static double complex r0; From 0fdc72b49d8f685fd7c6006eb6c026b3bdc0f7da Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Thu, 9 Jan 2025 18:30:16 +0900 Subject: [PATCH 16/17] improving comments --- LibHR/Observables/meson_measurements.c | 36 +++++++++++++++++++------- LibHR/Observables/polyakov.c | 13 +++++++--- LibHR/Observables/sources.c | 35 +++++++++++++++---------- 3 files changed, 58 insertions(+), 26 deletions(-) diff --git a/LibHR/Observables/meson_measurements.c b/LibHR/Observables/meson_measurements.c index ea9c20263..07959f7a4 100644 --- a/LibHR/Observables/meson_measurements.c +++ b/LibHR/Observables/meson_measurements.c @@ -768,8 +768,22 @@ void measure_formfactor_fixed(int ti, int tf, int dt, int nm, double *m, int n_m void create_smear_source_and_measure(spinor_field* source, spinor_field* prop, int t, int x, int y, int z, int nm, double* m, int n_mom, double epsilon_source, int Nsmear_source, int APE_N,int conf_num, char* label) { /** - * This function creates a smeared source and measures the meson correlators with point sink. - * The APE smearing on the gauge links are applied if `APE_N` is not 0. + * @brief Creates a Wuppetal smeared source and measures the meson correlators with point sink. + * The APE smearing on the gauge links are applied if `APE_N` is not 0. + * The smearing function, `smearing_function`, is imbedded in `create_smeared_source_with_APE` and `create_smeared_source`. + * See more details in that function. + * + * @param source The source field, $q(x)$. + * @param prop The fermion propagator for measuring the correlators. + * @param epsilon_source The smearing parameter, $\epsilon$, at the source. + * @param Nsmear_source The number of iterations for the Wuppertal smearing at the source. + * @param APE_N The number of iteration for the APE smearing, determining using `create_smeared_source_with_APE` or `create_smeared_source` + * @param t The time slice where the source is located. + * @param x The location of the source. + * @param y The location of the source. + * @param z The location of the source. + * @param nm Number of the measuring fermion mass + * @param m fermion masses */ for (int k = 0;k `pu_gauge_APE_f` + * @brief Applies the Wuppertal smearing to the source field with APE smeared gauge field, + * which is similar to the function `smearing_function` but changes the gauge field to the APE smeared one: + * `pu_gauge_f` -> `pu_gauge_APE_f` */ int ix, x, y, z, ix_up, ix_right, ix_front, ix_left, ix_back, ix_down; @@ -886,12 +886,17 @@ void smearing_function_with_APE(spinor_field *source, int tau, double epsilon){ void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear){ /** - * This function creates a smeared source field at a given position (t,x,y,z) with a given color by applying - * the Wuppertal smearing function based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. + * @brief Creates a smeared source field at a given position (t,x,y,z) with a given color by applying + * the Wuppertal smearing function on a point source based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. + * + * @param t The time slice where the source is located. + * @param x The location of the source. + * @param y The location of the source. + * @param z The location of the source. * @param Nsmear is the iteration number at the source - * @param epsilon is the step size. + * @param epsilon is the step size, $\epsilon$ in the reference. * - * The gauge field is the original one. + * @note The gauge field is the original one. */ int beta; @@ -913,13 +918,17 @@ void create_smeared_source(spinor_field *source, int t, int x, int y, int z, int void create_smeared_source_with_APE(spinor_field *source, int t, int x, int y, int z, int color, double epsilon, int Nsmear){ /** - * This function creates a smeared source field at a given position (t,x,y,z) with a given color by applying - * the Wuppertal smearing function based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. - * + * @brief Creates a smeared source field at a given position (t,x,y,z) with a given color by applying + * the Wuppertal smearing function on a point source based on Eq. (9) in https://arxiv.org/pdf/1602.05525.pdf. + * + * @param t The time slice where the source is located. + * @param x The location of the source. + * @param y The location of the source. + * @param z The location of the source. * @param Nsmear is the iteration number at the source - * @param epsilon is the step size. + * @param epsilon is the step size, $\epsilon$ in the reference. * - * The gauge field is the APE smeared one. + * @note The gauge field is the APE smeared one. */ int beta; From a2baa3f7b2574a260bfc2028c3c3234cb3f2387c Mon Sep 17 00:00:00 2001 From: Ho Hsiao <162536018+theHoHsiao@users.noreply.github.com> Date: Fri, 20 Jun 2025 15:20:29 +0900 Subject: [PATCH 17/17] add comments for `cooling_SPN` and helper functions --- LibHR/Utils/suN_utils.c | 68 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/LibHR/Utils/suN_utils.c b/LibHR/Utils/suN_utils.c index 9fc9cced3..61888f619 100644 --- a/LibHR/Utils/suN_utils.c +++ b/LibHR/Utils/suN_utils.c @@ -450,6 +450,20 @@ void covariant_project_to_suNg(suNg *u) void cooling_SPN(suNg* g_out, suNg* g_in, suNg* g_tilde, int cooling){ + /** + * @brief Cooling algorithm for SP(N) gauge fields applies the beta=infinity limit of the Cabibbo-Mariani algorithm, + * which gives the maximum value of Re Tr(g_tilde g_in). + * + * We follow the procedure described in https://arxiv.org/pdf/hep-lat/0404008, for the SU(N) case. + * + * @param g_out: Pointer to the output gauge field after cooling. + * @param g_in: Pointer to the input gauge field before cooling, which is U in Eq. (25) of the reference. + * @param g_tilde: Pointer to the input gauge field before cooling, which is \tilde{U} in Eq. (25) of the reference. + * @param cooling: Number of cooling iterations to perform. + * @note `subgrb` and `subgrb_tau` are helper functions to extract the SU(2) sub-matrices, + * while `vmxsu2` and `vmxsu2_tau` are helper functions to multiply the SU(2) sub-matrices. + * See detail explanations in Appendex A of the reference: https://arxiv.org/abs/2010.15781. + */ suNgfull B, U_tilde, S; _suNg_expand( U_tilde, *g_tilde); @@ -475,7 +489,20 @@ void cooling_SPN(suNg* g_out, suNg* g_in, suNg* g_tilde, int cooling){ } } + void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11){ + /** + * @brief Extracts an SU(2) subgroup from a Sp(N) matrix and applies a transformation. + * + * This function targets the SU(2) subgroup defined by the columns `i1col` and `i2col` + * of the input Sp(N) matrix `B11`. It constructs an SU(2) of type Eq. (A4) in the paper: https://arxiv.org/abs/2010.15781, + * normalizes it, and applies it to both `C11` and `B11` via left multiplication. + * + * @param i1col Index of the first column in Sp(N) matrix `B11`. + * @param i2col Index of the second column in Sp(N) matrix `B11`. + * @param B11 Pointer to the Sp(N) matrix to be updated in-place. + * @param C11 Pointer to a Sp(N) matrix also updated by the same SU(2) transformation. + */ double complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array double UMAG; // [ 0 1 ] @@ -509,7 +536,20 @@ void subgrb(int i1col, int i2col, suNgfull* B11, suNgfull* C11){ vmxsu2(i1col,i2col,B11,A11); } + void vmxsu2(int i1, int i2, suNgfull* A, double complex B[4]){ + + /** + * @brief Multiplies columns i1 and i2 of a Sp(N) matrix by a 2×2 SU(2) matrix. + * + * Applies an SU(2) rotation on columns i1 and i2 of matrix `A`, using matrix `B` + * provided in row-major order as a 1D array of 4 complex numbers. + * + * @param i1 First column index of Sp(N) matrix A. + * @param i2 Second column index of Sp(N) matrix A. + * @param A Pointer to the Sp(N) matrix. + * @param B SU(2) matrix in vector form: {B00, B01, B10, B11}. + */ double complex C[NG*2]; double complex ztmp1,ztmp2; @@ -530,7 +570,21 @@ void vmxsu2(int i1, int i2, suNgfull* A, double complex B[4]){ } } + void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11){ + + /** + * @brief Applies an SU(2) tau-subgroup transformation that mixes upper and lower blocks of Sp(N). + * + * Constructs a SU(2) matrix of type Eq. (A5) in the paper: https://arxiv.org/abs/2010.15781, + * + * Applies the resulting SU(2) matrix to both `B11` and `C11`. + * + * @param n1 Index in upper block. + * @param n2 Index in lower block. + * @param B11 Matrix B = S \tilde{U}^\dagger, used to generate the SU(2) transformation. + * @param C11 Matrix S, updated with the transformation. + */ double complex F11, F12, A11[4], ztmp1, ztmp2; // A11 is the extracted su2 matrix in 1-d array double UMAG; // [ 0 1 ] @@ -563,7 +617,21 @@ void subgrb_tau(int n1, int n2, suNgfull* B11, suNgfull* C11){ vmxsu2_tau(n1, n2, B11, A11); } + void vmxsu2_tau(int n1, int n2, suNgfull* A, complex B[4]){ + + /** + * @brief Applies a (A5)-type SU(2) transformation to a Sp(N) matrix. + * + * Performs the matrix update defined by a (A5)-type SU(2) matrix. + * A temporary copy is used to safely apply the transformation to four sub-blocks: + * (n1, n1), (n2, n2), (NG/2+n1, NG/2+n1), and (NG/2+n2, NG/2+n2). + * + * @param n1 Index in upper block. + * @param n2 Index in lower block. + * @param A Pointer to the Sp(N) matrix to be transformed. + * @param B SU(2) matrix in row-major form. + */ suNgfull C; double complex ztmp1, ztmp2;