Skip to content

Commit

Permalink
initialization of sources and propagator spinores + inversion + even-…
Browse files Browse the repository at this point in the history 
…odd lexic conversione
  • Loading branch information
@simone-romiti
simone-romiti committed Jun 16, 2023
1 parent f11e8da commit 20081dc
Showing 1 changed file with 132 additions and 52 deletions.
184 changes: 132 additions & 52 deletions meas/correlators.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,26 +260,15 @@ void light_correlators_measurement(const int traj, const int id, const int ieo)
*current), <S,S>, <S,P>, <P,P> for all the choices of h1 and h2 (see eq. 18 and 20 of
*https://arxiv.org/pdf/1005.2042.pdf)
*
* i1, i2 = source and sink operator indices from the list of the operators in the input file.
* i1, i2 = operator indices from the list of the operators in the input file.
*
*
******************************************************/
void heavy_quarks_correlators_measurement(const int traj, const int t0, const int ieo, const int i1,
const int i2) {
tm_stopwatch_push(&g_timers, __func__, ""); // timer for profiling
int i, j, t, tt, t0; // dummy indices
double *Cpp = NULL; // array of the values of the correlator C(t)
double res = 0.; // result of the accumulation of MPI partial sums
float tmp; // dummy variable
operator* optr1; // pointer to the operator (see some lines later)
#ifdef TM_USE_MPI
double mpi_res = 0., mpi_respa = 0., mpi_resp4 = 0.;
// send buffer for MPI_Gather
double *sCpp = NULL; //, *sCpa = NULL, *sCp4 = NULL;
#endif
FILE *ofs; // output file stream
char filename[TM_OMEAS_FILENAME_LENGTH]; // file path
spinor phi; // full spinor of 4 components (we use the spin dilution)

/* building the stichastic propagator spinors needed for the correlators */

init_operators(); // initialize the operators in the action (if not already done)
if (no_operators < 1) {
Expand All @@ -294,7 +283,8 @@ void heavy_quarks_correlators_measurement(const int traj, const int t0, const in
}

// selecting the operators i1 and i2 from the list of operators initialized before
optr1 = &operator_list[i1];
operator* optr1 = &operator_list[i1];
operator* optr2 = &operator_list[i2];

bool b1 = (optr1->type != DBTMWILSON && optr1->type != DBCLOVER);
if (b1) {
Expand Down Expand Up @@ -345,6 +335,129 @@ void heavy_quarks_correlators_measurement(const int traj, const int t0, const in
optr1->kappa, optr1->mubar, optr1->epsbar);
}


// even-odd spinor fields for the light and heavy doublet correlators
// source+propagator, even-odd, 2 flavors
// psi = psi[(s,p)][eo][f][x][alpha][c]
// Note: propagator in th sense that it is D^{-1}*source after the inversion
init_spinor_field(VOLUMEPLUSRAND / 2, 1); // initialize g_spinor_field so that I can copy it

This comment has been minimized.

Sign in to view

Copy link
@kostrzewa

kostrzewa Jun 16, 2023 •
edited
Loading

Member

@simone-romiti Be careful with this: this has already been called in hmc_tm or offline_measurement to initialize the global set of spinor field buffers.

For arrays of temporary spinor fields, you want to use allocate_spinor_field_array or init_solver_field.
spinor ****l_eo_spinor_field = (spinor ****)malloc(8 * VOLUME / 2);
spinor ****h_eo_spinor_field = (spinor ****)malloc(8 * VOLUME / 2);
// initalize to zero the source spinor fields (necessary??? isn't it enough to use the source
// generator?)
for (size_t i_sp = 0; i_sp < 2; i_sp++) { // source or propagator

for (size_t i_eo = 0; i_eo < 2; i_eo++) { // even-odd
for (size_t i_f = 0; i_f < 2; i++) { // up-down flavor
assign(l_eo_spinor_field[0][i_eo][i_f], g_spinor_field[0], VOLUME / 2);
assign(h_eo_spinor_field[0][i_eo][i_f], g_spinor_field[0], VOLUME / 2);
}
}
}

// initalize the random sources
// one source for each --> spin dilution
for (size_t i_f = 0; i_f < 2; i_f++) {
for (size_t src_d = 0; src_d < 4; src_d++) {
const unsigned int seed_i = src_d + measurement_list[id].seed;
// light doublet
eo_source_spinor_field_spin_diluted_oet_ts(l_eo_spinor_field[0][0][i_f],
l_eo_spinor_field[0][1][i_f], t0, src_d,
sample, traj, seed_i);
// heavy doublet
eo_source_spinor_field_spin_diluted_oet_ts(h_eo_spinor_field[0][0][i_f],
h_eo_spinor_field[0][1][i_f], t0, src_d,
sample, traj, seed_i);
}
}

// heavy doublet: for each even-odd block, I multiply by (1 + i*tau_2)
// the basis for the inversion should be the same as for the light
// --> I will multiply later by the inverse, namely at the end of the inversion
for (size_t i_eo = 0; i_eo < 2; i_eo++) {
for (size_t i_f = 0; i_f < 2; i_f++) {
// (u,d) --> [(1+i*tau_2)/sqrt(2)] * (u,d) , stored at temporarely in the propagator spinors (used as dummy spinors)
mul_one_pm_itau2(h_eo_spinor_field[1][i_eo][i_f], h_eo_spinor_field[1][i_eo][i_f+1],
h_eo_spinor_field[0][i_eo][i_f], h_eo_spinor_field[0][i_eo][i_f+1], +1.0, VOLUME / 2);
// assigning the result to the first components (the sources).
// The propagators will be overwritten with the inversion
assign(h_eo_spinor_field[0][i_eo][i_f], h_eo_spinor_field[1][i_eo][i_f], VOLUME / 2);
}
}

/*
assign the sources and propagator pointes for the operators
these need to be know when calling the inverter() method
*/

// light doublet
optr1->sr0 = l_eo_spinor_field[0][0][0];
optr1->sr1 = l_eo_spinor_field[0][0][1];
optr1->sr2 = l_eo_spinor_field[0][1][0];
optr1->sr3 = l_eo_spinor_field[0][1][1];

optr1->prop0 = l_eo_spinor_field[1][0][0];
optr1->prop1 = l_eo_spinor_field[1][0][1];
optr1->prop2 = l_eo_spinor_field[1][1][0];
optr1->prop3 = l_eo_spinor_field[1][1][1];

// heavy doublet
optr2->sr0 = h_eo_spinor_field[0][0][0];
optr2->sr1 = h_eo_spinor_field[0][0][1];
optr2->sr2 = h_eo_spinor_field[0][1][0];
optr2->sr3 = h_eo_spinor_field[0][1][1];

optr2->prop0 = h_eo_spinor_field[1][0][0];
optr2->prop1 = h_eo_spinor_field[1][0][1];
optr2->prop2 = h_eo_spinor_field[1][1][0];
optr2->prop3 = h_eo_spinor_field[1][1][1];

// inverting the Dirac operators for the light and heavy doublet
// op_id = i1 or i2, index_start = 0, write_prop = 0
optr1->inverter(i1, 0, 0);
optr2->inverter(i2, 0, 0);

// conclude the change of basis for the heavy doublet
for (size_t i_eo = 0; i_eo < 2; i_eo++) {
for (size_t i_f = 0; i_f < 2; i_f++) {
// (u,d) --> [(1+i*tau_2)/sqrt(2)] * (u,d) , stored at temporarely in the propagator spinors (used as dummy spinors)
mul_one_pm_itau2(h_eo_spinor_field[1][i_eo][i_f], h_eo_spinor_field[1][i_eo][i_f+1],
h_eo_spinor_field[0][i_eo][i_f], h_eo_spinor_field[0][i_eo][i_f+1], -1.0, VOLUME / 2);
// assigning the result to the first components (the sources).
// The propagators will be overwritten with the inversion
assign(h_eo_spinor_field[0][i_eo][i_f], h_eo_spinor_field[1][i_eo][i_f], VOLUME / 2);
}
}

// now we switch from even-odd representation to standard
// propagator, 2 flavors : psi = psi[f][x][alpha][c]
spinor ***l_propagator = (spinor ****)malloc(2 * VOLUME);
spinor ***h_propagator = (spinor ****)malloc(2 * VOLUME);
for (size_t i_f = 0; i_f < 2; i_f++) {
convert_eo_to_lexic(l_propagator[i_f], l_spinor_field[1][0][i_f], l_spinor_field[1][1][i_f]);
convert_eo_to_lexic(h_propagator[i_f], h_spinor_field[1][0][i_f], h_spinor_field[1][1][i_f]);
}

/*
Now that I have all the propagators (all in the basis of
https://arxiv.org/pdf/1005.2042.pdf) I can build the correlators of eq. 20
*/

spinor phi; // dummy spinor

// light-light

int i, j, t, tt, t0; // dummy indices
double *Cpp = NULL; // array of the values of the correlator C(t)
double res = 0.; // result of the accumulation of MPI partial sums
float tmp; // dummy variable
#ifdef TM_USE_MPI
double mpi_res = 0., mpi_respa = 0., mpi_resp4 = 0.;
// send buffer for MPI_Gather
double *sCpp = NULL; //, *sCpa = NULL, *sCp4 = NULL;
#endif
FILE *ofs; // output file stream
char filename[TM_OMEAS_FILENAME_LENGTH]; // file path
#ifdef TM_USE_MPI
sCpp = (double *)calloc(T, sizeof(double));
// sCpa = (double*) calloc(T, sizeof(double));
Expand All @@ -357,46 +470,13 @@ void heavy_quarks_correlators_measurement(const int traj, const int t0, const in
#else
Cpp = (double *)calloc(T, sizeof(double));
// Cpa = (double*) calloc(T, sizeof(double));
// Cp4 = (double*) calloc(T, sizeof(double));
#endif
// /// ??? what to use here
// source_generation_pion_only(g_spinor_field[0], g_spinor_field[1], t0, sample, traj,
// measurement_list[id].seed);
// initialize the random sources
for (size_t srd_d = 0; srd_d < 2; srd_d++) {
full_source_spinor_field_spin_diluted_oet_ts(g_bispinor_field[i1][src_d], t0, src_d, sample, traj,
measurement_list[id].seed);
}
// Cp4 = (double*) calloc(T, sizeof(double));// INIT SOME SPINOR FIELD FOR THE LIGHT

// ??? not sure how to use g_bispinor
spinor *up_spinors = &g_bispinor_field[i1][0][0];
spinor *down_spinors = &g_bispinor_field[0][1];
// sources
optr1->sr0 = g_bispinor_field[i1][0][0];
optr1->sr1 = g_bispinor_field[i1][0][1];
optr1->sr2 = down_spinors[0];
optr1->sr3 = down_spinors[1];
// propagators, i.e. D^{-1}*eta (eta = stochastic vector for the inversion)
optr1->prop0 = g_bispinor_field[i1][0][2];
optr1->prop1 = g_bispinor_field[i1][0][3];
optr1->prop2 = down_spinors[2];
optr1->prop3 = down_spinors[3];

// even-odd sites for the down
// sources
optr1->sr2 = g_bispinor_field[1][0];
optr1->sr3 = g_bispinor_field[1][1];
// propagators, i.e. D^{-1}*eta (eta = stochastic vector for the inversion)
optr1->prop2 = g_bispinor_field[1][2];
optr1->prop3 = g_bispinor_field[1][3];
#endif

// op_id = 0, index_start = 0, write_prop = 0
optr1->inverter(i1, 0, 0);
// heavy-light

/* now we bring it to normal format */
/* here we use implicitly DUM_MATRIX and DUM_MATRIX+1 */
convert_eo_to_lexic(down_spinors[DUM_MATRIX], down_spinors[2], down_spinors[3]);
convert_eo_to_lexic(down_spinors[DUM_MATRIX], down_spinors[2], down_spinors[3]);


/* now we sum only over local space for every t */
for (t = 0; t < T; t++) {
Expand Down

0 comments on commit 20081dc

Please sign in to comment.