From fd663253217bef8e83991cad8928d2a9a3056308 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 14 Dec 2018 17:39:11 +0000 Subject: [PATCH] pure QED test and copyright update --- Grid/qcd/action/gauge/Photon.h | 4 +- tests/core/Test_qed.cc | 138 +++++++++++++++++++++++++++++++++ 2 files changed, 141 insertions(+), 1 deletion(-) create mode 100644 tests/core/Test_qed.cc diff --git a/Grid/qcd/action/gauge/Photon.h b/Grid/qcd/action/gauge/Photon.h index e4a2fef4..f059fcf3 100644 --- a/Grid/qcd/action/gauge/Photon.h +++ b/Grid/qcd/action/gauge/Photon.h @@ -4,9 +4,11 @@ Source file: ./lib/qcd/action/gauge/Photon.h - Copyright (C) 2015 +Copyright (C) 2015-2018 Author: Peter Boyle + Author: Antonin Portelli + Author: James Harrison 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 diff --git a/tests/core/Test_qed.cc b/tests/core/Test_qed.cc new file mode 100644 index 00000000..f3b33454 --- /dev/null +++ b/tests/core/Test_qed.cc @@ -0,0 +1,138 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: tests/core/Test_qed.cc + +Copyright (C) 2015-2018 + +Author: Antonin Portelli +Author: James Harrison + +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 + +using namespace Grid; +using namespace QCD; + +typedef PeriodicGaugeImpl 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 + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(4, vComplex::Nsimd()); + std::vector 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 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 < 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.back(); ++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::avgWilsonLoop(expA, iw, iw) * 3; + logWlAvg[iw-1] -= 2*log(wlA); + wlA = WilsonLoops::avgTimelikeWilsonLoop(expA, iw, iw) * 3; + logWlTime[iw-1] -= 2*log(wlA); + wlA = WilsonLoops::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; +}