su2hmc
Loading...
Searching...
No Matches
Functions
integrate.c File Reference

Symplectic integrator code. More...

#include <su2hmc.h>
Include dependency graph for integrate.c:

Go to the source code of this file.

Functions

int Gauge_Update (const double d, double *pp, Complex *ut[2], Complex_f *ut_f[2])
 Gauge update for the integration step of the HMC.
 
int Momentum_Update (const double d, const double *dSdpi, double *pp)
 Wrapper for the momentum update during the integration step of the HMC.
 
int Leapfrog (Complex *ut[2], Complex_f *ut_f[2], Complex *X0, Complex *X1, Complex *Phi, double *dk[2], float *dk_f[2], double *dSdpi, double *pp, int *iu, int *id, Complex *gamval, Complex_f *gamval_f, int *gamin, Complex jqq, float beta, float akappa, int stepl, float dt, double *ancg, int *itot, float proby)
 Leapfrog integrator. Each trajectory step takes the form of p->p+dt/2,u->u+dt,p->p+dt/2 In practice this is implemented for the entire trajectory as p->p+dt/2,u->u+dt,p->p+dt,u->u+dt,p->p+dt,...p->p+dt/2,u->u+dt,p->p+dt/2.
 
int OMF2 (Complex *ut[2], Complex_f *ut_f[2], Complex *X0, Complex *X1, Complex *Phi, double *dk[2], float *dk_f[2], double *dSdpi, double *pp, int *iu, int *id, Complex *gamval, Complex_f *gamval_f, int *gamin, Complex jqq, float beta, float akappa, int stepl, float dt, double *ancg, int *itot, float proby)
 OMF second order five step integrator.
 
int OMF4 (Complex *ut[2], Complex_f *ut_f[2], Complex *X0, Complex *X1, Complex *Phi, double *dk[2], float *dk_f[2], double *dSdpi, double *pp, int *iu, int *id, Complex *gamval, Complex_f *gamval_f, int *gamin, Complex jqq, float beta, float akappa, int stepl, float dt, double *ancg, int *itot, float proby)
 OMF fourth order eleven step integrator.
 

Detailed Description

Symplectic integrator code.

Definition in file integrate.c.

Function Documentation

◆ Gauge_Update()

int Gauge_Update ( const double d,
double * pp,
Complex * ut[2],
Complex_f * ut_f[2] )

Gauge update for the integration step of the HMC.

Parameters
dtGauge step size
ppMomentum field
utDouble precision gauge fields
ut_fSingle precision gauge fields
Returns
Zero on success, integer error code otherwise.

Definition at line 7 of file integrate.c.

7 {
8 /*
9 * @brief Generates new trial fields
10 *
11 * @see cuGauge_Update (CUDA Wrapper)
12 *
13 * @param d: Half lattice spacing
14 * @param pp: Momentum field
15 * @param ut[0]: First colour field
16 * @param ut[1]: Second colour field
17 *
18 * @returns Zero on success, integer error code otherwise
19 */
20 char *funcname = "Gauge_Update";
21#ifdef __NVCC__
22 cuGauge_Update(d,pp,ut[0],ut[1],dimGrid,dimBlock);
23#else
24#pragma omp parallel for simd collapse(2) aligned(pp:AVX)
25 for(int i=0;i<kvol;i++)
26 for(int mu = 0; mu<ndim; mu++){
27 /*
28 * Sticking to what was in the FORTRAN for variable names.
29 * CCC for cosine SSS for sine AAA for...
30 * Re-exponentiating the force field. Can be done analytically in SU(2)
31 * using sine and cosine which is nice
32 */
33
34 double AAA = d*sqrt(pp[i*nadj*ndim+mu]*pp[i*nadj*ndim+mu]\
35 +pp[(i*nadj+1)*ndim+mu]*pp[(i*nadj+1)*ndim+mu]\
36 +pp[(i*nadj+2)*ndim+mu]*pp[(i*nadj+2)*ndim+mu]);
37 double CCC = cos(AAA);
38 double SSS = d*sin(AAA)/AAA;
39 Complex a11 = CCC+I*SSS*pp[(i*nadj+2)*ndim+mu];
40 Complex a12 = pp[(i*nadj+1)*ndim+mu]*SSS + I*SSS*pp[i*nadj*ndim+mu];
41 //b11 and b12 are ut[0] and ut[1] terms, so we'll use ut[1] directly
42 //but use b11 for ut[0] to prevent RAW dependency
43 complex b11 = ut[0][i*ndim+mu];
44 ut[0][i*ndim+mu] = a11*b11-a12*conj(ut[1][i*ndim+mu]);
45 ut[1][i*ndim+mu] = a11*ut[1][i*ndim+mu]+a12*conj(b11);
46 }
47#endif
48 Reunitarise(ut);
49 //Get trial fields from accelerator for halo exchange
50 Trial_Exchange(ut,ut_f);
51 return 0;
52}
Trial_Exchange
int Trial_Exchange(Complex *ut[2], Complex_f *ut_f[2])
Exchanges the trial fields.
Definition par_mpi.c:1178
nadj
#define nadj
adjacent spatial indices
Definition sizes.h:175
kvol
#define kvol
Sublattice volume.
Definition sizes.h:154
Complex
#define Complex
Double precision complex number.
Definition sizes.h:58
ndim
#define ndim
Dimensions.
Definition sizes.h:179
Reunitarise
int Reunitarise(Complex *ut[2])
Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1.
Definition matrices.c:1012

References Complex, Complex_f, kvol, nadj, ndim, Reunitarise(), and Trial_Exchange().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Leapfrog()

int Leapfrog ( Complex * ut[2],
Complex_f * ut_f[2],
Complex * X0,
Complex * X1,
Complex * Phi,
double * dk[2],
float * dk_f[2],
double * dSdpi,
double * pp,
int * iu,
int * id,
Complex * gamval,
Complex_f * gamval_f,
int * gamin,
Complex jqq,
float beta,
float akappa,
int stepl,
float dt,
double * ancg,
int * itot,
float proby )

Leapfrog integrator. Each trajectory step takes the form of p->p+dt/2,u->u+dt,p->p+dt/2 In practice this is implemented for the entire trajectory as p->p+dt/2,u->u+dt,p->p+dt,u->u+dt,p->p+dt,...p->p+dt/2,u->u+dt,p->p+dt/2.

Parameters
utDouble precision colour fields
ut_fSingle precision colour fields
X0Up/down partitioned pseudofermion field
X1Holder for the partitioned fermion field, then the conjugate gradient output
PhiPseudofermion field
dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
dk_f\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
dSdpiThe force
ppMomentum field
iu,idLattice indices
gaminGamma indices
gamvalDouble precision gamma matrices rescaled by kappa
gamval_fSingle precision gamma matrices rescaled by kappa
jqqDiquark source
akappaHopping parameter
betaInverse gauge coupling
steplSteps per trajectory
dtStep size
ancgCounter for average conjugate gradient iterations
itotTotal average conjugate gradient iterations
probyTermination probability for random trajectory length
Returns
Zero on success, integer error code otherwise

Definition at line 65 of file integrate.c.

68{
69 //This was originally in the half-step of the FORTRAN code, but it makes more sense to declare
70 //it outside the loop. Since it's always being subtracted we'll define it as negative
71 const double d =-dt*0.5;
72 //Half step forward for p
73 //=======================
74#ifdef _DEBUG
75 printf("Evaluating force on rank %i\n", rank);
76#endif
77 Force(dSdpi, 1, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
78 Momentum_Update(d,dSdpi,pp);
79 //Main loop for classical time evolution
80 //======================================
81 bool end_traj=false; int step =1;
82 // for(int step = 1; step<=stepmax; step++){
83 do{
84#ifdef _DEBUG
85 if(!rank)
86 printf("Step: %d\n", step);
87#endif
88 //The FORTRAN redefines d=dt here, which makes sense if you have a limited line length.
89 //I'll stick to using dt though.
90 //step (i) st(t+dt)=st(t)+p(t+dt/2)*dt;
91 //Note that we are moving from kernel to kernel within the default streams so don't need a Device_Sync here
92 Gauge_Update(dt,pp,ut,ut_f);
93
94 //p(t+3et/2)=p(t+dt/2)-dSds(t+dt)*dt
95 // Force(dSdpi, 0, rescgg);
96 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
97
98 // if(step>=stepl*4.0/5.0 && (step>=stepl*(6.0/5.0) || Par_granf()<proby)){
99 if(step==stepl){
100 //Final trajectory has a half momentum step
101 Momentum_Update(d,dSdpi,pp);
102 *itot+=step;
103 *ancg/=step;
104 end_traj=true;
105 break;
106 }
107 else{
108 //Otherwise, there's a half step at the end and start of each trajectory, so we combine them into one full step.
109 Momentum_Update(-dt,dSdpi,pp);
110 step++;
111 }
112 }while(!end_traj);
113
114 return 0;
115}
Gauge_Update
int Gauge_Update(const double d, double *pp, Complex *ut[2], Complex_f *ut_f[2])
Gauge update for the integration step of the HMC.
Definition integrate.c:7
Momentum_Update
int Momentum_Update(const double d, const double *dSdpi, double *pp)
Wrapper for the momentum update during the integration step of the HMC.
Definition integrate.c:53
rank
int rank
The MPI rank.
Definition par_mpi.c:22
rescgg
#define rescgg
Conjugate gradient residue for update.
Definition sizes.h:240
Force
int Force(double *dSdpi, int iflag, double res1, Complex *X0, Complex *X1, Complex *Phi, Complex *ut[2], Complex_f *ut_f[2], unsigned int *iu, unsigned int *id, Complex *gamval, Complex_f *gamval_f, int *gamin, double *dk[2], float *dk_f[2], Complex_f jqq, float akappa, float beta, double *ancg)
Calculates the force at each intermediate time.
Definition force.c:94

References Complex, Complex_f, Force(), Gauge_Update(), Momentum_Update(), rank, and rescgg.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ Momentum_Update()

int Momentum_Update ( const double d,
const double * dSdpi,
double * pp )
inline

Wrapper for the momentum update during the integration step of the HMC.

Parameters
dtGauge step size
ppMomentum field
dSdpiForce field
Returns
Zero on success, integer error code otherwise.

Definition at line 53 of file integrate.c.

54{
55#ifdef __NVCC__
56 cublasDaxpy(cublas_handle,kmom, &d, dSdpi, 1, pp, 1);
57#elif defined USE_BLAS
58 cblas_daxpy(kmom, d, dSdpi, 1, pp, 1);
59#else
60#pragma omp parallel for simd
61 for(int i=0;i<kmom;i++)
62 pp[i]+=d*dSdpi[i];
63#endif
64}
kmom
#define kmom
sublattice momentum sites
Definition sizes.h:184

References kmom.

Here is the caller graph for this function:

◆ OMF2()

int OMF2 ( Complex * ut[2],
Complex_f * ut_f[2],
Complex * X0,
Complex * X1,
Complex * Phi,
double * dk[2],
float * dk_f[2],
double * dSdpi,
double * pp,
int * iu,
int * id,
Complex * gamval,
Complex_f * gamval_f,
int * gamin,
Complex jqq,
float beta,
float akappa,
int stepl,
float dt,
double * ancg,
int * itot,
float proby )

OMF second order five step integrator.

Parameters
utDouble precision colour fields
ut_fSingle precision colour fields
X0Up/down partitioned pseudofermion field
X1Holder for the partitioned fermion field, then the conjugate gradient output
PhiPseudofermion field
dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
dk_f\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
dSdpiThe force
ppMomentum field
iu,idLattice indices
gaminGamma indices
gamvalDouble precision gamma matrices rescaled by kappa
gamval_fSingle precision gamma matrices rescaled by kappa
jqqDiquark source
akappaHopping parameter
betaInverse gauge coupling
steplSteps per trajectory
dtStep size
ancgCounter for average conjugate gradient iterations
itotTotal average conjugate gradient iterations
probyTermination probability for random trajectory length
Returns
Zero on success, integer error code otherwise

Definition at line 116 of file integrate.c.

119{
120 const double lambda=0.5-(pow(2.0*sqrt(326.0)+36.0,1.0/3.0)/12.0)+1.0/(6*pow(2.0*sqrt(326.0) + 36.0,1.0/3.0));
121 //const double lambda=1.0/6.0;
122 // const double lambda=0.5;
123
124 //Gauge update by half dt
125 const double dU = dt*0.5;
126
127 //Momentum updates by lambda, 2lambda and (1-2lambda) in the middle
128 const double dp= -lambda*dt;
129 const double dp2= 2.0*dp;
130 const double dpm= -(1.0-2.0*lambda)*dt;
131 //Initial step forward for p
132 //=======================
133#ifdef _DEBUG
134 printf("Evaluating force on rank %i\n", rank);
135#endif
136 Force(dSdpi, 1, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
137#ifdef _DEBUG
138 printf("Initial force dSdpi[0]=%e\n", dSdpi[0]);
139#endif
140 //Initial momentum update
141 Momentum_Update(dp,dSdpi,pp);
142#ifdef _DEBUG
143 printf("Initial OMF2 momentum pp[0]=%e\n", pp[0]);
144#endif
145
146 //Main loop for classical time evolution
147 //======================================
148 bool end_traj=false; int step =1;
149 // for(int step = 1; step<=stepmax; step++){
150 do{
151#ifdef _DEBUG
152 if(!rank)
153 printf("Step: %d\n", step);
154#endif
155 //First gauge update
156 Gauge_Update(dU,pp,ut,ut_f);
157
158 //Calculate force for middle momentum update
159 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
160 //Now do the middle momentum update
161 Momentum_Update(dpm,dSdpi,pp);
162#ifdef _DEBUG
163 printf("Middle OMF2 momentum pp[0]=%e\n", pp[0]);
164 #endif
165
166 //Second gauge update
167 Gauge_Update(dU,pp,ut,ut_f);
168
169 //Calculate force for second momentum update
170 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
171
172 //if(step>=stepl*4.0/5.0 && (step>=stepl*(6.0/5.0) || Par_granf()<proby)){
173 if(step==stepl){
174 //Final momentum step
175 Momentum_Update(dp,dSdpi,pp);
176#ifdef _DEBUG
177 printf("Final OMF2 momentum pp[0]=%e\n", pp[0]);
178 #endif
179 *itot+=step;
180 //Two force terms, so an extra factor of two in the average?
181 //Or leave it as it was, to get the average CG iterations per trajectory rather than force
182 *ancg/=step;
183 end_traj=true;
184 break;
185 }
186 else{
187 //Since we apply the momentum at the start and end of a step we instead apply a double step here
188 Momentum_Update(dp2,dSdpi,pp);
189 step++;
190 }
191 }while(!end_traj);
192 return 0;
193}

References Complex, Complex_f, Force(), Gauge_Update(), Momentum_Update(), rank, and rescgg.

Here is the call graph for this function:
Here is the caller graph for this function:

◆ OMF4()

int OMF4 ( Complex * ut[2],
Complex_f * ut_f[2],
Complex * X0,
Complex * X1,
Complex * Phi,
double * dk[2],
float * dk_f[2],
double * dSdpi,
double * pp,
int * iu,
int * id,
Complex * gamval,
Complex_f * gamval_f,
int * gamin,
Complex jqq,
float beta,
float akappa,
int stepl,
float dt,
double * ancg,
int * itot,
float proby )

OMF fourth order eleven step integrator.

Parameters
utDouble precision colour fields
ut_fSingle precision colour fields
X0Up/down partitioned pseudofermion field
X1Holder for the partitioned fermion field, then the conjugate gradient output
PhiPseudofermion field
dk\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\)
dk_f\(\left(1+\gamma_0\right)e^{-\mu}\) and \(\left(1-\gamma_0\right)e^\mu\) float
dSdpiThe force
ppMomentum field
iu,idLattice indices
gaminGamma indices
gamvalDouble precision gamma matrices rescaled by kappa
gamval_fSingle precision gamma matrices rescaled by kappa
jqqDiquark source
akappaHopping parameter
betaInverse gauge coupling
steplSteps per trajectory
dtStep size
ancgCounter for average conjugate gradient iterations
itotTotal average conjugate gradient iterations
probyTermination probability for random trajectory length
Returns
Zero on success, integer error code otherwise

Momentum updates

Outer updates depend on r1. We have two of these, doubled for between full steps

Middle updates depend on r3

Inner updates depend on r1 and r3

Gauge updates. These depend on r2 and r4

Outer gauge update depends on r2

Middle gauge update depends on r4

Inner gauge update depends on r2 and r4

Definition at line 194 of file integrate.c.

197{
198 //These values were lifted from openqcd-fastsum, and should probably be tuned for QC2D. They also probably never
199 //will be...
200 const double r1 = 0.08398315262876693;
201 const double r2 = 0.2539785108410595;
202 const double r3 = 0.6822365335719091;
203 const double r4 = -0.03230286765269967;
206 const double dpO= -r1*dt;
207 const double dpO2= 2*dpO;
209 const double dpM= -r3*dt;
211 const double dpI= -(0.5-r1-r3)*dt;
212
215 const double duO = dt*r2;
217 const double duM = dt*r4;
219 const double duI = dt*(1-2*(r2+r4));
220
221 //Initial step forward for p
222 //=======================
223#ifdef _DEBUG
224 printf("Evaluating force on rank %i\n", rank);
225#endif
226 Force(dSdpi, 1, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
227 Momentum_Update(dpO,dSdpi,pp);
228
229 //Main loop for classical time evolution
230 //======================================
231 bool end_traj=false; int step =1;
232 // for(int step = 1; step<=stepmax; step++){
233 do{
234#ifdef _DEBUG
235 if(!rank)
236 printf("Step: %d\n", step);
237#endif
238 //First outer gauge update
239 Gauge_Update(duO,pp,ut,ut_f);
240
241 //Calculate force for first middle momentum update
242 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
243 //Now do the first middle momentum update
244 Momentum_Update(dpM,dSdpi,pp);
245
246 //First middle gauge update
247 Gauge_Update(duM,pp,ut,ut_f);
248
249 //Calculate force for first inner momentum update
250 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
251 //Now do the first inner momentum update
252 Momentum_Update(dpI,dSdpi,pp);
253
254 //Inner gauge update
255 Gauge_Update(duI,pp,ut,ut_f);
256
257 //Calculate force for second inner momentum update
258 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
259 //Now do the second inner momentum update
260 Momentum_Update(dpI,dSdpi,pp);
261
262 //Second middle gauge update
263 Gauge_Update(duM,pp,ut,ut_f);
264
265 //Calculate force for second middle momentum update
266 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
267 //Now do the second middle momentum update
268 Momentum_Update(dpM,dSdpi,pp);
269
270 //Second outer gauge update
271 Gauge_Update(duO,pp,ut,ut_f);
272
273 //Calculate force for outer momentum update
274 Force(dSdpi, 0, rescgg,X0,X1,Phi,ut,ut_f,iu,id,gamval,gamval_f,gamin,dk,dk_f,jqq,akappa,beta,ancg);
275
276 //Outer momentum update depends on if we've finished the trajectory
277 //if(step>=stepl*4.0/5.0 && (step>=stepl*(6.0/5.0) || Par_granf()<proby)){
278 if(step==stepl){
279 //Final momentum step
280 Momentum_Update(dpO,dSdpi,pp);
281 *itot+=step;
282
283 //Four force terms, so an extra factor of four in the average?
284 //Or leave it as it was, to get the average CG iterations per trajectory rather than force
285 *ancg/=step;
286 end_traj=true;
287 break;
288 }
289 else{
290 //Since we apply the momentum at the start and end of a step we instead apply a double step here
291 Momentum_Update(dpO2,dSdpi,pp);
292 step++;
293 }
294 }while(!end_traj);
295 return 0;
296}

References Complex, Complex_f, Force(), Gauge_Update(), Momentum_Update(), rank, and rescgg.

Here is the call graph for this function:
Here is the caller graph for this function: