1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-14 09:45:36 +00:00
Grid/tests/core/Test_qed.cc

138 lines
4.8 KiB
C++

/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/core/Test_qed.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: James Harrison <J.Harrison@soton.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
#include <Grid/Grid.h>
using namespace Grid;
typedef PeriodicGaugeImpl<QedGImplR> QedPeriodicGImplR;
typedef PhotonR::GaugeField EmField;
typedef PhotonR::GaugeLinkField EmComp;
const int NCONFIGS = 20;
const int NWILSON = 10;
int main(int argc, char *argv[])
{
// initialization
Grid_init(&argc, &argv);
std::cout << GridLogMessage << "Grid initialized" << std::endl;
// QED stuff
Coordinate latt_size = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(4, vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
GridCartesian grid(latt_size,simd_layout,mpi_layout);
GridParallelRNG pRNG(&grid);
PhotonR photon(&grid, PhotonR::Gauge::coulomb, PhotonR::ZmScheme::qedL);
EmField a(&grid);
EmField expA(&grid);
Complex imag_unit(0, 1);
Real wlA;
std::vector<Real> logWlAvg(NWILSON, 0.0), logWlTime(NWILSON, 0.0), logWlSpace(NWILSON, 0.0);
pRNG.SeedFixedIntegers({1, 2, 3, 4});
std::cout << GridLogMessage << "Wilson loop calculation beginning" << std::endl;
for(int ic = 0; ic < NCONFIGS; ic++){
std::cout << GridLogMessage << "Configuration " << ic <<std::endl;
photon.StochasticField(a, pRNG);
// Exponentiate photon field
expA = exp(imag_unit*a);
// Calculate zero-modes
std::vector<EmField::vector_object::scalar_object> zm;
std::cout << GridLogMessage << "Total zero-mode norm 2 "
<< std::sqrt(norm2(sum(a))) << std::endl;
std::cout << GridLogMessage << "Spatial zero-mode norm 2" << std::endl;
sliceSum(a, zm, grid.Nd() - 1);
for (unsigned int t = 0; t < latt_size.size(); ++t)
{
std::cout << GridLogMessage << "t = " << t << " " << std::sqrt(norm2(zm[t])) << std::endl;
}
// Calculate divergence
EmComp diva(&grid), amu(&grid);
diva = Zero();
for (unsigned int mu = 0; mu < grid.Nd(); ++mu)
{
amu = peekLorentz(a, mu);
diva += amu - Cshift(amu, mu, -1);
if (mu == grid.Nd() - 2)
{
std::cout << GridLogMessage << "Spatial divergence norm 2 " << std::sqrt(norm2(diva)) << std::endl;
}
}
std::cout << GridLogMessage << "Total divergence norm 2 " << std::sqrt(norm2(diva)) << std::endl;
// Calculate Wilson loops
for(int iw=1; iw<=NWILSON; iw++){
wlA = WilsonLoops<QedPeriodicGImplR>::avgWilsonLoop(expA, iw, iw) * 3;
logWlAvg[iw-1] -= 2*log(wlA);
wlA = WilsonLoops<QedPeriodicGImplR>::avgTimelikeWilsonLoop(expA, iw, iw) * 3;
logWlTime[iw-1] -= 2*log(wlA);
wlA = WilsonLoops<QedPeriodicGImplR>::avgSpatialWilsonLoop(expA, iw, iw) * 3;
logWlSpace[iw-1] -= 2*log(wlA);
}
}
std::cout << GridLogMessage << "Wilson loop calculation completed" << std::endl;
// Calculate Wilson loops
// From A. Portelli's PhD thesis:
// size -2*log(W)
// 1 0.500000000(1)
// 2 1.369311535(1)
// 3 2.305193057(1)
// 4 3.261483854(1)
// 5 4.228829967(1)
// 6 5.203604529(1)
// 7 6.183728249(1)
// 8 7.167859805(1)
// 9 8.155091868(1)
// 10 9.144788116(1)
for(int iw=1; iw<=10; iw++){
std::cout << GridLogMessage << iw << 'x' << iw << " Wilson loop" << std::endl;
std::cout << GridLogMessage << "-2*log(W) average: " << logWlAvg[iw-1]/NCONFIGS << std::endl;
std::cout << GridLogMessage << "-2*log(W) timelike: " << logWlTime[iw-1]/NCONFIGS << std::endl;
std::cout << GridLogMessage << "-2*log(W) spatial: " << logWlSpace[iw-1]/NCONFIGS << std::endl;
}
// epilogue
std::cout << GridLogMessage << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}