diff --git a/tests/forces/Test_dwf_force_eofa.cc b/tests/forces/Test_dwf_force_eofa.cc index 80d36934..525178d0 100644 --- a/tests/forces/Test_dwf_force_eofa.cc +++ b/tests/forces/Test_dwf_force_eofa.cc @@ -86,7 +86,9 @@ int main (int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); + RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_dwf_gpforce.cc b/tests/forces/Test_dwf_gpforce.cc index 28133cc6..1fa1c6e4 100644 --- a/tests/forces/Test_dwf_gpforce.cc +++ b/tests/forces/Test_dwf_gpforce.cc @@ -84,6 +84,13 @@ int main (int argc, char ** argv) GparityDomainWallFermionR::ImplParams params; params.twists = twists; + /* + params.boundary_phases[0] = 1.0; + params.boundary_phases[1] = 1.0; + params.boundary_phases[2] = 1.0; + params.boundary_phases[3] =- 1.0; + */ + GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); Dw.M (phi,Mphi); @@ -96,6 +103,16 @@ int main (int argc, char ** argv) Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + + // ***************************************************************************************** + // *** There is a funny negative sign in all derivatives. This is - UdSdU. *** + // *** *** + // *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign *** + // *** UdSdU is negated relative to what I think - call what is returned mUdSdU, *** + // *** and insert minus sign *** + // ***************************************************************************************** + + UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy. FermionField Ftmp (FGrid); @@ -106,7 +123,7 @@ int main (int argc, char ** argv) RealD Hmom = 0.0; RealD Hmomprime = 0.0; LatticeColourMatrix mommu(UGrid); - LatticeColourMatrix forcemu(UGrid); + LatticeColourMatrix mUdSdUmu(UGrid); LatticeGaugeField mom(UGrid); LatticeGaugeField Uprime(UGrid); @@ -114,10 +131,20 @@ int main (int argc, char ** argv) SU::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg - Hmom -= real(sum(trace(mommu*mommu))); + // Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR + // + // Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu); + // Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra + // P is i P^a_\mu T^a, not Pa Ta + // + // Integrator.h: H = Hmom + sum S(action) + Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR; PokeIndex(mom,mommu,mu); + // -- Drops factor of "i" in the U update: U' = e^{P dt} U [ _not_ e^{iPdt}U ]. P is anti hermitian already + // -- Udot = p U + // fourth order exponential approx autoView( mom_v, mom, CpuRead); autoView( U_v , U, CpuRead); @@ -134,8 +161,8 @@ int main (int argc, char ** argv) ; }); } - std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <(mom,mu); - std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); - std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); - mommu=Ta(mommu)*2.0; + mommu=Ta(mommu); // projectForce , GaugeImplTypes.h PokeIndex(UdSdU,mommu,mu); } for(int mu=0;mu(mom,mu); - std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<(UdSdU,mu); - std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<(UdSdU,mu); + mUdSdUmu= PeekIndex(UdSdU,mu); mommu = PeekIndex(mom,mu); - // Update PF action density - dS = dS+trace(mommu*forcemu)*dt; + // + // Derive HMC eom: + // + // Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0 + // + // + // Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0 + // + // EOM: + // + // pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above + // + // dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit + // + // dH_mom/dt = - 2 trace (p pdot)/Denom + // + // dH_tot / dt = 0 <= pdot = - Denom * mUdSdU + // + // dH_mom/dt = 2 trace (p mUdSdU ) + // + // True Momentum delta H has a dt^2 piece + // + // dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom + // = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f). + // = dSmom + dSmom2 + // - dSmom = dSmom - trace(mommu*forcemu) * dt; - dSmom2 = dSmom2 - trace(forcemu*forcemu) *(0.25* dt*dt); + dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x. - // Update mom action density - mommu = mommu + forcemu*(dt*0.5); + dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2 + + dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant - Hmomprime -= real(sum(trace(mommu*mommu))); + // Update mom action density . Verbatim update_P in Integrator.h + mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;; + + Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR; } @@ -199,20 +233,25 @@ int main (int argc, char ** argv) ComplexD dSm = sum(dSmom); ComplexD dSm2 = sum(dSmom2); + std::cout << GridLogMessage <<"dSm "<< dSm< CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_laplacian_force.cc b/tests/forces/Test_laplacian_force.cc index 18508860..dbaf1cbd 100644 --- a/tests/forces/Test_laplacian_force.cc +++ b/tests/forces/Test_laplacian_force.cc @@ -46,6 +46,7 @@ int main (int argc, char ** argv) std::vector seeds({1,2,3,4}); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers({4,5,6,7}); GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(std::vector({15,91,21,3})); @@ -67,7 +68,7 @@ int main (int argc, char ** argv) LaplacianAdjointField Laplacian(&Grid, CG, LapPar, Kappa); GeneralisedMomenta LaplacianMomenta(&Grid, Laplacian); LaplacianMomenta.M.ImportGauge(U); - LaplacianMomenta.MomentaDistribution(pRNG);// fills the Momenta with the correct distr + LaplacianMomenta.MomentaDistribution(sRNG,pRNG);// fills the Momenta with the correct distr std::cout << std::setprecision(15); diff --git a/tests/forces/Test_mobius_force.cc b/tests/forces/Test_mobius_force.cc index ba7bc363..d2326a81 100644 --- a/tests/forces/Test_mobius_force.cc +++ b/tests/forces/Test_mobius_force.cc @@ -69,7 +69,14 @@ int main (int argc, char ** argv) RealD M5=1.8; RealD b=0.5; RealD c=0.5; - MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + + WilsonImplParams p; + p.boundary_phases[0] = 1.0; + p.boundary_phases[1] = 1.0; + p.boundary_phases[2] = 1.0; + p.boundary_phases[3] =- 1.0; + + MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,p); Ddwf.M (phi,Mphi); ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p @@ -82,24 +89,44 @@ int main (int argc, char ** argv) Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp); + // ***************************************************************************************** + // *** There is a funny negative sign in all derivatives. This is - UdSdU. *** + // *** *** + // *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign *** + // *** UdSdU is negated relative to what I think - call what is returned mUdSdU, *** + // *** and insert minus sign *** + // ***************************************************************************************** + + UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy. + LatticeFermion Ftmp (FGrid); //////////////////////////////////// // Modify the gauge field a little //////////////////////////////////// - RealD dt = 0.0001; + RealD dt = 0.001; + RealD Hmom = 0.0; + RealD Hmomprime = 0.0; LatticeColourMatrix mommu(UGrid); - LatticeColourMatrix forcemu(UGrid); + LatticeColourMatrix mUdSdUmu(UGrid); LatticeGaugeField mom(UGrid); LatticeGaugeField Uprime(UGrid); for(int mu=0;mu::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg - PokeIndex(mom,mommu,mu); + // Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR + // + // Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu); + // Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra + // P is i P^a_\mu T^a, not Pa Ta + // + // Integrator.h: H = Hmom + sum S(action) + Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR; + // fourth order exponential approx autoView( U_v , U, CpuRead); autoView( mom_v, mom, CpuRead); @@ -115,6 +142,7 @@ int main (int argc, char ** argv) ; }); } + std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <(UdSdU,mu); - mommu=Ta(mommu)*2.0; + mommu=Ta(mommu); PokeIndex(UdSdU,mommu,mu); } for(int mu=0;mu(UdSdU,mu); + + mUdSdUmu= PeekIndex(UdSdU,mu); mommu = PeekIndex(mom,mu); - // Update PF action density - dS = dS+trace(mommu*forcemu)*dt; + // + // Derive HMC eom: + // + // Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0 + // + // + // Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0 + // + // EOM: + // + // pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above + // + // dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit + // + // dH_mom/dt = - 2 trace (p pdot)/Denom + // + // dH_tot / dt = 0 <= pdot = - Denom * mUdSdU + // + // dH_mom/dt = 2 trace (p mUdSdU ) + // + // True Momentum delta H has a dt^2 piece + // + // dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom + // = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f). + // = dSmom + dSmom2 + // + + dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x. + + dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2 + + dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant + + mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;; + + Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR; } ComplexD dSpred = sum(dS); + ComplexD dSm = sum(dSmom); + ComplexD dSm2 = sum(dSmom2); - std::cout << GridLogMessage << " -- S "< CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U" diff --git a/tests/forces/Test_mobius_gpforce_eofa.cc b/tests/forces/Test_mobius_gpforce_eofa.cc index 9c80b2aa..7f114615 100644 --- a/tests/forces/Test_mobius_gpforce_eofa.cc +++ b/tests/forces/Test_mobius_gpforce_eofa.cc @@ -93,7 +93,8 @@ int main (int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false); - Meofa.refresh(U, RNG5); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4); + Meofa.refresh(U, sRNG, RNG5 ); RealD S = Meofa.S(U); // pdag M p // get the deriv of phidag M phi with respect to "U"