diff --git a/Include/Utils/gaussian_smearing.h b/Include/Utils/gaussian_smearing.h new file mode 100644 index 000000000..5223db826 --- /dev/null +++ b/Include/Utils/gaussian_smearing.h @@ -0,0 +1,16 @@ +#ifndef GAUSS_SMEARING +#define GAUSS_SMEARING + +#include "spinor_field.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void gaussian_smearing(spinor_field *restrict out, spinor_field *restrict in, suNf_field *gauge_f, double alpha); + + +#ifdef __cplusplus +} +#endif +#endif \ No newline at end of file diff --git a/Include/utils.h b/Include/utils.h index 0e987010d..4ff743cef 100644 --- a/Include/utils.h +++ b/Include/utils.h @@ -13,6 +13,7 @@ #include "Utils/eva_deflation.h" #include "Utils/gauge_anisotropy.h" #include "Utils/gaugefix.h" +#include "Utils/gaussian_smearing.h" #include "Utils/HYP_smearing.h" #include "Utils/mat_utils.h" #include "Utils/eigval.h" @@ -31,4 +32,4 @@ // void cross_prod(suNg_vector *v1, suNg_vector *v2, suNg_vector *v3); //TOOD: this is not defined in libhr // void cross_prod_flt(suNg_vector_flt *v1, suNg_vector_flt *v2, suNg_vector_flt *v3); //TOOD: this is not defined in libhr -#endif +#endif \ No newline at end of file diff --git a/LibHR/Utils/gaussian_smearing.c b/LibHR/Utils/gaussian_smearing.c new file mode 100644 index 000000000..2511968f0 --- /dev/null +++ b/LibHR/Utils/gaussian_smearing.c @@ -0,0 +1,105 @@ +/****************************************************************************** +* File: gaussian_smearing.c +* +* Gaussian smearing for fermions +* +* Author: Pietro Butti +* Started: Nov. 2025 +******************************************************************************/ + +#include "utils.h" +#include "libhr_core.h" +#include "memory.h" +#include "IO/logger.h" +#include "inverters.h" +#include "update.h" + + +/*******************************************\ + out <= F * in = (1 + aH) * in + + H(n,m)psi(m) = sum_i [U_i(n) psi(n+i) + U_i(n-i)^* psi(n-i)] +\*******************************************/ +void gaussian_smearing(spinor_field *restrict out, spinor_field *restrict in, suNf_field *gauge_f, double alpha) { +#ifdef CHECK_SPINOR_MATCHING + error((in == NULL) || (out == NULL), 1, "gaussian_smearing [gaussian_smearing.c]", "Attempt to access unallocated memory space"); + error(in == out, 1, "gaussian_smearing [gaussian_smearing.c]", "Input and output fields must be different"); + error(out->type == &glat_even && in->type == &glat_even, 1, "gaussian_smearing [gaussian_smearing.c]", "Spinors don't match! (1)"); + error(out->type == &glat_odd && in->type == &glat_odd, 1, "gaussian_smearing [gaussian_smearing.c]", "Spinors don't match! (2)"); +#endif + + /************************ loop over all lattice sites *************************/ + /* start communication of input spinor field */ + start_sendrecv_spinor_field(in); + + _PIECE_FOR(out->type, ixp) { +#ifdef WITH_MPI + if (ixp == out->type->inner_master_pieces) { + /* wait for spinor to be transfered */ + complete_sendrecv_spinor_field(in); + } +#endif + + _SITE_FOR(out->type, ixp, ix) { +#if defined(_OPENMP) && defined(WITH_PROBE_MPI) + register int thread0 = hr_threadId(); +#endif + + int i_dw, i_up; + double factor = 1. / (1. + 6. * alpha); + suNf *u, *u_dw; + suNf_spinor *psi_up, *psi_dw, *r, *rold; + suNf_vector chi0, chi1, chi2, chi3; + + r = _FIELD_AT(out, ix); + rold = _FIELD_AT(in, ix); + + _spinor_zero_f(*r); + + for (int idir = 1; idir < 4; ++idir) { + i_dw = idn(ix, idir); + i_up = iup(ix, idir); + + u = (gauge_f->ptr) + coord_to_index(ix, idir); + u_dw = (gauge_f->ptr) + coord_to_index(i_dw, idir); + + psi_up = _FIELD_AT(in, i_up); + psi_dw = _FIELD_AT(in, i_dw); + + // up ------------------------------------------ + // chi_i <= U_ij(n) * psi_j(n+dir) + _suNf_multiply(chi0, *(u), (*psi_up).c[0]); + _suNf_multiply(chi1, *(u), (*psi_up).c[1]); + _suNf_multiply(chi2, *(u), (*psi_up).c[2]); + _suNf_multiply(chi3, *(u), (*psi_up).c[3]); + + // r <= r + chi_i + _vector_add_assign_f((*r).c[0], chi0); + _vector_add_assign_f((*r).c[1], chi1); + _vector_add_assign_f((*r).c[2], chi2); + _vector_add_assign_f((*r).c[3], chi3); + + // dw ------------------------------------------ + // chi_i <= U_ij^dag(n-dir) * psi_j(n-dir) + _suNf_inverse_multiply(chi0, *(u_dw), (*psi_dw).c[0]); + _suNf_inverse_multiply(chi1, *(u_dw), (*psi_dw).c[1]); + _suNf_inverse_multiply(chi2, *(u_dw), (*psi_dw).c[2]); + _suNf_inverse_multiply(chi3, *(u_dw), (*psi_dw).c[3]); + + // r <= r + chi_i + _vector_add_assign_f((*r).c[0], chi0); + _vector_add_assign_f((*r).c[1], chi1); + _vector_add_assign_f((*r).c[2], chi2); + _vector_add_assign_f((*r).c[3], chi3); + } + + /******************************** end of loop *********************************/ + _spinor_lc_f(*r, factor, *rold, alpha * factor, *r); + +#ifdef WITH_PROBE_MPI + if (thread0 == 0) { probe_mpi(); } +#endif + + } /* SITE_FOR */ + } /* PIECE FOR */ +} \ No newline at end of file diff --git a/TestProgram/Sources/check_smeared_stoch_sources.c b/TestProgram/Sources/check_smeared_stoch_sources.c new file mode 100644 index 000000000..f0505ab54 --- /dev/null +++ b/TestProgram/Sources/check_smeared_stoch_sources.c @@ -0,0 +1,280 @@ +/****************************************************************************** +* File: check_smeared_stoch_sources.c +* +* Checks of smeared stochastic sources +* +* Author: Pietro Butti +* Started: Sept. 2025 +******************************************************************************/ + +#include "libhr.h" + +/* Mesons parameters */ +typedef struct input_mesons { + char mstring[256]; + + /* for the reading function */ + input_record_t read[5]; + double csw; + double alpha; //PIETRO// + int nhits_2pt; + +} input_mesons; + +#define init_input_mesons(varname) \ + { \ + .read = { \ + { "fermion mass", "mes:mass = %s", STRING_T, (varname).mstring }, \ + { "csw coefficient", "mes:csw = %lf", DOUBLE_T, &(varname).csw }, \ + { "smearing alpha", "mes:alpha = %lf", DOUBLE_T, &(varname).alpha }, \ + { "number of noisy sources per cnfg for 2pt fn", "mes:nhits_2pt = %d", INT_T, &(varname).nhits_2pt }, \ + { NULL, NULL, INT_T, NULL } \ + } \ + } + +input_glb glb_ip = init_input_glb(glb_ip); +input_mesons mes_ip = init_input_mesons(mes_ip); + +void create_alltimes_pt_source_rand(spinor_field *source, int *locs) { + int t, beta, col, lidx, ix; + double ran; + + int Xsrc, Ysrc, Zsrc; + + for (beta = 0; beta < 4; ++beta) { +#ifdef WITH_GPU + zero_spinor_field_cpu(&source[beta]); +#endif + zero_spinor_field(&source[beta]); + } + + if (PID == 0) { + // Generate random spatial location + for (t = 0; t < GLB_T; t++) { + locs[4 * t] = t; + ranlxd(&ran, 1); + locs[1 + 4 * t] = (int)(ran * GLB_X); + ranlxd(&ran, 1); + locs[2 + 4 * t] = (int)(ran * GLB_Y); + ranlxd(&ran, 1); + locs[3 + 4 * t] = (int)(ran * GLB_Z); + } + } // end rank 0 global task +#ifdef WITH_MPI + // Broadcast info to every processes + MPI_Bcast(locs, 4 * GLB_T, MPI_INT, 0, GLB_COMM); +#endif + + for (t = 0; t < T; t++) { + Xsrc = locs[1 + 4 * (zerocoord[0] + t)]; + Ysrc = locs[2 + 4 * (zerocoord[0] + t)]; + Zsrc = locs[3 + 4 * (zerocoord[0] + t)]; + + if ((zerocoord[1] <= Xsrc && Xsrc < zerocoord[1] + X) && (zerocoord[2] <= Ysrc && Ysrc < zerocoord[2] + Y) && + (zerocoord[3] <= Zsrc && Zsrc < zerocoord[3] + Z)) { + // Select point + ix = ipt(t, Xsrc - zerocoord[1], Ysrc - zerocoord[2], Zsrc - zerocoord[3]); + + // Fill source + for (col = 0; col < NF; col++) { + for (beta = 0; beta < 4; beta++) { + lidx = beta + col * 4; + _FIELD_AT(&source[lidx], ix)->c[beta].c[col] = 1.; + } + } + } + } + + // ----------- + + for (beta = 0; beta < 4; ++beta) { +#ifdef WITH_GPU + copy_to_gpu(source + beta); +#endif + start_sendrecv_spinor_field(source + beta); + complete_sendrecv_spinor_field(source + beta); + } +} + +void create_alltimes_pt_source_loc(spinor_field *source, int Xsrc, int Ysrc, int Zsrc) { + int t, beta, col, lidx, ix; + + for (beta = 0; beta < 4; ++beta) { +#ifdef WITH_GPU + zero_spinor_field_cpu(&source[beta]); +#endif + zero_spinor_field(&source[beta]); + } + + // Check whether (t,xs,ys,zs) is in this process + if (zerocoord[1] <= Xsrc && Xsrc < zerocoord[1] + X && zerocoord[2] <= Ysrc && Ysrc < zerocoord[2] + Y && + zerocoord[3] <= Zsrc && Zsrc < zerocoord[3] + Z) { + for (t = 0; t < T; t++) { + // Select point + ix = ipt(t, Xsrc - zerocoord[1], Ysrc - zerocoord[2], Zsrc - zerocoord[3]); + + // Fill source + for (col = 0; col < NF; col++) { + for (beta = 0; beta < 4; beta++) { + lidx = beta + col * 4; + _FIELD_AT(&source[lidx], ix)->c[beta].c[col] = 1.; + } + } + } + } + + for (beta = 0; beta < 4; ++beta) { +#ifdef WITH_GPU + copy_to_gpu(source + beta); +#endif + start_sendrecv_spinor_field(source + beta); + complete_sendrecv_spinor_field(source + beta); + } +} + +double smear_factor(int vecp[3], double alpha) { + double cosp = cos(2. * M_PI * (double)vecp[0] / (double)GLB_X) + cos(2. * M_PI * (double)vecp[1] / (double)GLB_Y) + + cos(2. * M_PI * (double)vecp[2] / (double)GLB_Z); + double z = 1. / (1. + 6. * alpha) * (1. + 2. * alpha * cosp); + + return z; +} + +hr_complex exp_i_p_x(int p0, int p1, int p2, int x0, int x1, int x2) { + return (cos(2. * M_PI / (double)GLB_X * (double)(p0 * x0)) + I * sin(2. * M_PI / (double)GLB_X * (double)(p0 * x0))) * + (cos(2. * M_PI / (double)GLB_Y * (double)(p1 * x1)) + I * sin(2. * M_PI / (double)GLB_Y * (double)(p1 * x1))) * + (cos(2. * M_PI / (double)GLB_Z * (double)(p2 * x2)) + I * sin(2. * M_PI / (double)GLB_Z * (double)(p2 * x2))); +} + +int check_alltime_rand(spinor_field *source, hr_complex *acc, int *vecp, int *locs, double alpha) { + int t, ix, iy, iz, ii, col, beta; + + hr_complex pw, tmp, accumulator; + // hr_complex f_p = smear_factor(vecp, alpha); + spinor_field *lsource; + + for (t = 0; t < T; t++) { + // Fourier transfrom + spin/color sum in a block + accumulator = 0. + I * 0.; + for (ix = 0; ix < X; ix++) { + for (iy = 0; iy < Y; iy++) { + for (iz = 0; iz < Z; iz++) { + // Select point + ii = ipt(t, ix, iy, iz); + + // Compute plane waves + pw = exp_i_p_x(-vecp[0], -vecp[1], -vecp[2], ix + zerocoord[1] - locs[1 + 4 * (zerocoord[0] + t)], + iy + zerocoord[2] - locs[2 + 4 * (zerocoord[0] + t)], + iz + zerocoord[3] - locs[3 + 4 * (zerocoord[0] + t)]); + + // Cycle over spin/color components + for (int ls = 0; ls < NF * 4; ls++) { + lsource = source + ls; + + for (col = 0; col < NF; col++) { + for (beta = 0; beta < 4; beta++) { + // int lidx = beta + col * 4; + // Accumulate + tmp = pw * _FIELD_AT(lsource, ii)->c[beta].c[col]; + // if (creal(tmp) * creal(tmp) + cimag(tmp) * cimag(tmp) > 1.e-12) { + // lprintf("not zero", 0, "point %d %d %d %d %d %d value %.10e %.10e\n", ix, iy, iz, t, col, + // beta, creal(tmp), cimag(tmp)); + // } + + accumulator += tmp; + } + } + } + } + } + } + // Accumulate + acc[t + zerocoord[0]] = accumulator; + } + + global_sum((double *)acc, 2 * GLB_T); + return 0; +} + +int main(int argc, char *argv[]) { + int n_sources; + int beta, col, lidx; + + hr_complex *acc; + + logger_map("DEBUG", "debug"); + logger_setlevel(0, 200); + + /* setup process id and communications */ + setup_process(&argc, &argv); + // read_input(pars_ip.read, get_input_filename()); + read_input(mes_ip.read, get_input_filename()); + + setup_gauge_fields(); + unit_u(u_gauge); + represent_gauge_field(); + + lprintf("CORR", 0, "smearing parameter : alpha = %e \n", mes_ip.alpha); + + n_sources = NF * 4; + // spinor_field *source = alloc_spinor_field(n_sources, &glattice); + // spinor_field *source_pt = alloc_spinor_field(n_sources, &glattice); + + // // Create one source per spin and color index + // create_alltimes_pt_source_loc(source_pt, xs, ys, zs); + + int *src_loc = malloc(4 * GLB_T * sizeof(int)); + spinor_field *src = alloc_spinor_field(n_sources, &glattice); + spinor_field *src_pt = alloc_spinor_field(n_sources, &glattice); + + create_alltimes_pt_source_rand(src_pt, src_loc); + for (int t = 0; t < GLB_T; t++) { + lprintf("MAIN", 0, "t=%i [%i, %i, %i]\n", src_loc[4 * t], src_loc[1 + 4 * t], src_loc[2 + 4 * t], src_loc[3 + 4 * t]); + } + + // Smear the source + for (col = 0; col < NF; col++) { + for (beta = 0; beta < 4; beta++) { + lidx = beta + col * 4; + // Fphi_cpu_(&source[lidx], &source_pt[lidx], mes_ip.alpha); + gaussian_smearing(&src[lidx], &src_pt[lidx], u_gauge_f, mes_ip.alpha); + } + } + + // Perform this to a representative of momenta ----------------------- + int moms[3 * 4] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1 }; + hr_complex f_p; + + for (int imom = 0; imom < 4; imom++) { + int vecp[3] = { moms[3 * imom], moms[1 + 3 * imom], moms[2 + 3 * imom] }; + f_p = smear_factor(vecp, mes_ip.alpha) * (double)n_sources; + + lprintf("TEST", 0, "========= p = [%i, %i, %i], f(p) = %e + i %e ========= \n", vecp[0], vecp[1], vecp[2], creal(f_p), + cimag(f_p)); + + // Run test + acc = (hr_complex *)malloc(GLB_T * sizeof(hr_complex)); + for (int t=0; t