diff --git a/.gitignore b/.gitignore index 40156f9d..94e866e2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ +# Doxygen stuff +html/* +latex/* + # Compiled Object files # ######################### *.slo diff --git a/Grid/log/Log.h b/Grid/log/Log.h index 2d663a3c..370b0428 100644 --- a/Grid/log/Log.h +++ b/Grid/log/Log.h @@ -179,11 +179,11 @@ extern GridLogger GridLogSolver; extern GridLogger GridLogError; extern GridLogger GridLogWarning; extern GridLogger GridLogMessage; -extern GridLogger GridLogDebug ; +extern GridLogger GridLogDebug; extern GridLogger GridLogPerformance; extern GridLogger GridLogDslash; -extern GridLogger GridLogIterative ; -extern GridLogger GridLogIntegrator ; +extern GridLogger GridLogIterative; +extern GridLogger GridLogIntegrator; extern GridLogger GridLogHMC; extern GridLogger GridLogMemory; extern GridLogger GridLogTracing; @@ -191,6 +191,41 @@ extern Colours GridLogColours; std::string demangle(const char* name) ; +template +inline std::string sjoin(Args&&... args) noexcept { + std::ostringstream msg; + (msg << ... << args); + return msg.str(); +} + +/*! @brief make log messages work like python print */ +template +inline void Grid_log(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << GridLogMessage << msg << std::endl; +} + +/*! @brief make warning messages work like python print */ +template +inline void Grid_warn(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << "\033[33m" << GridLogWarning << msg << "\033[0m" << std::endl; +} + +/*! @brief make error messages work like python print */ +template +inline void Grid_error(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << "\033[31m" << GridLogError << msg << "\033[0m" << std::endl; +} + +/*! @brief make pass messages work like python print */ +template +inline void Grid_pass(Args&&... args) { + std::string msg = sjoin(std::forward(args)...); + std::cout << "\033[32m" << GridLogMessage << msg << "\033[0m" << std::endl; +} + #define _NBACKTRACE (256) extern void * Grid_backtrace_buffer[_NBACKTRACE]; diff --git a/Grid/qcd/smearing/HISQSmearing.h b/Grid/qcd/smearing/HISQSmearing.h new file mode 100644 index 00000000..529ea090 --- /dev/null +++ b/Grid/qcd/smearing/HISQSmearing.h @@ -0,0 +1,389 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/smearing/HISQSmearing.h + +Copyright (C) 2023 + +Author: D. A. Clarke + +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 +*************************************************************************************/ +/* + @file HISQSmearing.h + @brief Declares classes related to HISQ smearing +*/ + + +#pragma once +#include +#include +#include + + +NAMESPACE_BEGIN(Grid); + + +// TODO: find a way to fold this into the stencil header. need to access grid to get +// Nd, since you don't want to inherit from QCD.h +/*! @brief append arbitrary shift path to shifts */ +template +void appendShift(std::vector& shifts, int dir, Args... args) { + Coordinate shift(Nd,0); + generalShift(shift, dir, args...); + // push_back creates an element at the end of shifts and + // assigns the data in the argument to it. + shifts.push_back(shift); +} + + +/*! @brief figure out the stencil index from mu and nu */ +accelerator_inline int stencilIndex(int mu, int nu) { + // Nshifts depends on how you built the stencil + int Nshifts = 6; + return Nshifts*nu + Nd*Nshifts*mu; +} + + +/*! @brief structure holding the link treatment */ +struct SmearingParameters{ + SmearingParameters(){} + Real c_1; // 1 link + Real c_naik; // Naik term + Real c_3; // 3 link + Real c_5; // 5 link + Real c_7; // 7 link + Real c_lp; // 5 link Lepage + SmearingParameters(Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) + : c_1(c1), + c_naik(cnaik), + c_3(c3), + c_5(c5), + c_7(c7), + c_lp(clp){} +}; + + +/*! @brief create fat links from link variables */ +template +class Smear_HISQ : public Gimpl { + +private: + GridCartesian* const _grid; + SmearingParameters _linkTreatment; + +public: + + INHERIT_GIMPL_TYPES(Gimpl); + typedef typename Gimpl::GaugeField GF; + typedef typename Gimpl::GaugeLinkField LF; + typedef typename Gimpl::ComplexField CF; + + // Don't allow default values here. + Smear_HISQ(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) + : _grid(grid), + _linkTreatment(c1,cnaik,c3,c5,c7,clp) { + assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); + assert(Nd == 4 && "HISQ smearing only defined for Nd==4"); + } + + // Allow to pass a pointer to a C-style, double array for MILC convenience + Smear_HISQ(GridCartesian* grid, double* coeff) + : _grid(grid), + _linkTreatment(coeff[0],coeff[1],coeff[2],coeff[3],coeff[4],coeff[5]) { + assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3"); + assert(Nd == 4 && "HISQ smearing only defined for Nd==4"); + } + + ~Smear_HISQ() {} + + // Intent: OUT--u_smr, u_naik + // IN--u_thin + void smear(GF& u_smr, GF& u_naik, GF& u_thin) const { + + SmearingParameters lt = this->_linkTreatment; + auto grid = this->_grid; + + // Create a padded cell of extra padding depth=1 and fill the padding. + int depth = 1; + PaddedCell Ghost(depth,grid); + GF Ughost = Ghost.Exchange(u_thin); + + // This is where auxiliary N-link fields and the final smear will be stored. + GF Ughost_fat(Ughost.Grid()); + GF Ughost_3link(Ughost.Grid()); + GF Ughost_5linkA(Ughost.Grid()); + GF Ughost_5linkB(Ughost.Grid()); + + // mu-nu plane stencil. We allow mu==nu to make indexing the stencil easier, + // but these entries will not be used. + std::vector shifts; + for(int mu=0;mu_permute,Nd)) U3matrix; + + int Nsites = U_v.size(); + auto gStencil_v = gStencil.View(); + + accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 3-link constructs + stencilElement SE0, SE1, SE2, SE3, SE4, SE5; + U3matrix U0, U1, U2, U3, U4, U5, W; + for(int nu=0;nu_offset; + SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu = SE1->_offset; + SE2 = gStencil_v.GetEntry(s+2,site); int x = SE2->_offset; + SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset; + SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu = SE4->_offset; + SE5 = gStencil_v.GetEntry(s+5,site); int x_m_mu = SE5->_offset; + + // When you're deciding whether to take an adjoint, the question is: how is the + // stored link oriented compared to the one you want? If I imagine myself travelling + // with the to-be-updated link, I have two possible, alternative 3-link paths I can + // take, one starting by going to the left, the other starting by going to the right. + U0 = coalescedReadGeneralPermute(U_v[x_p_mu ](nu),SE0->_permute,Nd); + U1 = coalescedReadGeneralPermute(U_v[x_p_nu ](mu),SE1->_permute,Nd); + U2 = coalescedReadGeneralPermute(U_v[x ](nu),SE2->_permute,Nd); + U3 = coalescedReadGeneralPermute(U_v[x_p_mu_m_nu](nu),SE3->_permute,Nd); + U4 = coalescedReadGeneralPermute(U_v[x_m_nu ](mu),SE4->_permute,Nd); + U5 = coalescedReadGeneralPermute(U_v[x_m_nu ](nu),SE4->_permute,Nd); + + // "left" "right" + W = U2*U1*adj(U0) + adj(U5)*U4*U3; + + // Save 3-link construct for later and add to smeared field. + coalescedWrite(U_3link_v[x](nu), W); + + // The index operator (x) returns the coalesced read on GPU. The view [] index returns + // a reference to the vector object. The [x](mu) returns a reference to the densely + // packed (contiguous in memory) mu-th element of the vector object. On CPU, + // coalescedRead/Write is the identity mapping assigning vector object to vector object. + // But on GPU it's non-trivial and maps scalar object to vector object and vice versa. + coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_3*W); + } + }) + + accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 5-link + stencilElement SE0, SE1, SE2, SE3, SE4, SE5; + U3matrix U0, U1, U2, U3, U4, U5, W; + int sigmaIndex = 0; + for(int nu=0;nu_offset; + SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu = SE1->_offset; + SE2 = gStencil_v.GetEntry(s+2,site); int x = SE2->_offset; + SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset; + SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu = SE4->_offset; + + U0 = coalescedReadGeneralPermute( U_v[x_p_mu ](nu ),SE0->_permute,Nd); + U1 = coalescedReadGeneralPermute(U_3link_v[x_p_nu ](rho),SE1->_permute,Nd); + U2 = coalescedReadGeneralPermute( U_v[x ](nu ),SE2->_permute,Nd); + U3 = coalescedReadGeneralPermute( U_v[x_p_mu_m_nu](nu ),SE3->_permute,Nd); + U4 = coalescedReadGeneralPermute(U_3link_v[x_m_nu ](rho),SE4->_permute,Nd); + U5 = coalescedReadGeneralPermute( U_v[x_m_nu ](nu ),SE4->_permute,Nd); + + W = U2*U1*adj(U0) + adj(U5)*U4*U3; + + if(sigmaIndex<3) { + coalescedWrite(U_5linkA_v[x](rho), W); + } else { + coalescedWrite(U_5linkB_v[x](rho), W); + } + + coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_5*W); + sigmaIndex++; + } + } + }) + + accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 7-link + stencilElement SE0, SE1, SE2, SE3, SE4, SE5; + U3matrix U0, U1, U2, U3, U4, U5, W; + int sigmaIndex = 0; + for(int nu=0;nu_offset; + SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu = SE1->_offset; + SE2 = gStencil_v.GetEntry(s+2,site); int x = SE2->_offset; + SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset; + SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu = SE4->_offset; + + U0 = coalescedReadGeneralPermute(U_v[x_p_mu](nu),SE0->_permute,Nd); + if(sigmaIndex<3) { + U1 = coalescedReadGeneralPermute(U_5linkB_v[x_p_nu](rho),SE1->_permute,Nd); + } else { + U1 = coalescedReadGeneralPermute(U_5linkA_v[x_p_nu](rho),SE1->_permute,Nd); + } + U2 = coalescedReadGeneralPermute(U_v[x](nu),SE2->_permute,Nd); + U3 = coalescedReadGeneralPermute(U_v[x_p_mu_m_nu](nu),SE3->_permute,Nd); + if(sigmaIndex<3) { + U4 = coalescedReadGeneralPermute(U_5linkB_v[x_m_nu](rho),SE4->_permute,Nd); + } else { + U4 = coalescedReadGeneralPermute(U_5linkA_v[x_m_nu](rho),SE4->_permute,Nd); + } + U5 = coalescedReadGeneralPermute(U_v[x_m_nu](nu),SE4->_permute,Nd); + + W = U2*U1*adj(U0) + adj(U5)*U4*U3; + + coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_7*W); + sigmaIndex++; + } + } + }) + + } // end mu loop + + // c1, c3, c5, c7 construct contributions + u_smr = Ghost.Extract(Ughost_fat) + lt.c_1*u_thin; + + // Load up U and V std::vectors to access thin and smeared links. + std::vector U(Nd, grid); + std::vector V(Nd, grid); + std::vector Vnaik(Nd, grid); + for (int mu = 0; mu < Nd; mu++) { + U[mu] = PeekIndex(u_thin, mu); + V[mu] = PeekIndex(u_smr, mu); + } + + for(int mu=0;mu(u_smr , V[mu] , mu); + PokeIndex(u_naik, Vnaik[mu], mu); + } + }; + + + // Intent: OUT--u_proj + // IN--u_mu + void projectU3(GF& u_proj, GF& u_mu) const { + + auto grid = this->_grid; + + LF V(grid), Q(grid), sqrtQinv(grid), id_3(grid), diff(grid); + CF c0(grid), c1(grid), c2(grid), g0(grid), g1(grid), g2(grid), S(grid), R(grid), theta(grid), + u(grid), v(grid), w(grid), den(grid), f0(grid), f1(grid), f2(grid); + + // Follow MILC 10.1103/PhysRevD.82.074501, eqs (B2-B3) and (C1-C8) + for (int mu = 0; mu < Nd; mu++) { + V = PeekIndex(u_mu, mu); + Q = adj(V)*V; + c0 = real(trace(Q)); + c1 = (1/2.)*real(trace(Q*Q)); + c2 = (1/3.)*real(trace(Q*Q*Q)); + S = (1/3.)*c1-(1/18.)*c0*c0; + if (norm2(S)<1e-28) { + g0 = (1/3.)*c0; g1 = g0; g2 = g1; + } else { + R = (1/2.)*c2-(1/3. )*c0*c1+(1/27.)*c0*c0*c0; + theta = acos(R*pow(S,-1.5)); + g0 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta-2*M_PI/3.); + g1 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta ); + g2 = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta+2*M_PI/3.); + } +// if (fabs(Q.determinant()/(g0*g1*g2)-1.0) > 1e-5) { SVD } + u = sqrt(g0) + sqrt(g1) + sqrt(g2); + v = sqrt(g0*g1) + sqrt(g0*g2) + sqrt(g1*g2); + w = sqrt(g0*g1*g2); + den = w*(u*v-w); + f0 = (-w*(u*u+v)+u*v*v)/den; + f1 = (-w-u*u*u+2.*u*v)/den; + f2 = u/den; + id_3 = 1.; + + sqrtQinv = f0*id_3 + f1*Q + f2*Q*Q; + + PokeIndex(u_proj, V*sqrtQinv, mu); + } + }; + + +// void derivative(const GaugeField& Gauge) const { +// }; +}; + + +NAMESPACE_END(Grid); \ No newline at end of file diff --git a/Grid/qcd/smearing/Smearing.h b/Grid/qcd/smearing/Smearing.h index da5ede72..41a305ae 100644 --- a/Grid/qcd/smearing/Smearing.h +++ b/Grid/qcd/smearing/Smearing.h @@ -5,4 +5,5 @@ #include #include #include +#include diff --git a/Grid/stencil/GeneralLocalStencil.h b/Grid/stencil/GeneralLocalStencil.h index 90fff953..9af9b834 100644 --- a/Grid/stencil/GeneralLocalStencil.h +++ b/Grid/stencil/GeneralLocalStencil.h @@ -137,5 +137,55 @@ public: }; + +//////////////////////////////////////////////// +// Some machinery to streamline making a stencil +//////////////////////////////////////////////// + +class shiftSignal { +public: + enum { + BACKWARD_CONST = 16, + NO_SHIFT = -1 + }; +}; + +// TODO: put a check somewhere that BACKWARD_CONST > Nd! + +/*! @brief signals that you want to go backwards in direction dir */ +inline int Back(const int dir) { + // generalShift will use BACKWARD_CONST to determine whether we step forward or + // backward. Trick inspired by SIMULATeQCD. + return dir + shiftSignal::BACKWARD_CONST; +} + +/*! @brief shift one unit in direction dir */ +template +void generalShift(Coordinate& shift, int dir) { + if (dir >= shiftSignal::BACKWARD_CONST) { + dir -= shiftSignal::BACKWARD_CONST; + shift[dir]+=-1; + } else if (dir == shiftSignal::NO_SHIFT) { + ; // do nothing + } else { + shift[dir]+=1; + } +} + +/*! @brief follow a path of directions, shifting one unit in each direction */ +template +void generalShift(Coordinate& shift, int dir, Args... args) { + if (dir >= shiftSignal::BACKWARD_CONST) { + dir -= shiftSignal::BACKWARD_CONST; + shift[dir]+=-1; + } else if (dir == shiftSignal::NO_SHIFT) { + ; // do nothing + } else { + shift[dir]+=1; + } + generalShift(shift, args...); +} + + NAMESPACE_END(Grid); diff --git a/examples/Example_plaquette.cc b/examples/Example_plaquette.cc new file mode 100644 index 00000000..faf17d82 --- /dev/null +++ b/examples/Example_plaquette.cc @@ -0,0 +1,183 @@ +/* + * Example_plaquette.cc + * + * D. Clarke + * + * Here I just want to create an incredibly simple main to get started with GRID and get used + * to its syntax. If the reader is like me, they vaguely understand something about lattice coding, + * they don't know a ton of C++, don't know much of the fine details, and certainly know nothing about GRID. + * + * Once you've made a new executable, like this one, you can bootstrap.sh again. At this point, + * the code should be able to find your new executable. You can tell that bootstrap.sh worked by + * having a look at Make.inc. You should see your executable inside there. + * + * Warning: This code illustrative only, not well tested, and not meant for production use. The best + * way to read this code is to start at the main. + * + */ + + +// All your mains should have this +#include +using namespace Grid; + + +// This copies what already exists in WilsonLoops.h. The point here is to be pedagogical and explain in +// detail what everything does so we can see how GRID works. +template class WLoops : public Gimpl { +public: + // Gimpl seems to be an arbitrary class. Within this class, it is expected that certain types are + // already defined, things like Scalar and Field. This macro includes a bunch of #typedefs that + // implement this equivalence at compile time. + INHERIT_GIMPL_TYPES(Gimpl); + + // Some example Gimpls can be found in GaugeImplementations.h, at the bottom. These are in turn built + // out of GaugeImplTypes, which can be found in GaugeImplTypes.h. The GaugeImplTypes contain the base + // field/vector/link/whatever types. These inherit from iScalar, iVector, and iMatrix objects, which + // are sort of the building blocks for gerenal math objects. The "i" at the beginning of these names + // indicates that they should be for internal use only. It seems like these base types have the + // acceleration, e.g. SIMD or GPU or what-have-you, abstracted away. How you accelerate these things + // appears to be controlled through a template parameter called vtype. + + // The general math/physics objects, such as a color matrix, are built up by nesting these objects. + // For instance a general color matrix has two color indices, so it's built up like + // iScalar &U, const int mu, const int nu) { + // Calls like CovShiftForward and CovShiftBackward have 3 arguments, and they multiply together + // the first and last argument. (Second arg gives the shift direction.) The CovShiftIdentityBackward + // has meanwhile only two arguments; it just returns the shifted (adjoint since backward) link. + plaq = Gimpl::CovShiftForward(U[mu],mu, + // Means Link*Cshift(field,mu,1), arguments are Link, mu, field in that order. + Gimpl::CovShiftForward(U[nu],nu, + Gimpl::CovShiftBackward(U[mu],mu, + // This means Cshift(adj(Link), mu, -1) + Gimpl::CovShiftIdentityBackward(U[nu], nu)))); + } + + // tr U_mu_nu(x) + static void traceDirPlaquette(ComplexField &plaq, const std::vector &U, const int mu, const int nu) { + // This .Grid() syntax seems to get the pointer to the GridBase. Apparently this is needed as argument + // to instantiate a Lattice object. + GaugeMat sp(U[0].Grid()); + dirPlaquette(sp, U, mu, nu); + plaq = trace(sp); + } + + // sum_mu_nu tr U_mu_nu(x) + static void sitePlaquette(ComplexField &Plaq, const std::vector &U) { + ComplexField sitePlaq(U[0].Grid()); + Plaq = Zero(); + // Nd=4 and Nc=3 are set as global constants in QCD.h + for (int mu = 1; mu < Nd; mu++) { + for (int nu = 0; nu < mu; nu++) { + traceDirPlaquette(sitePlaq, U, mu, nu); + Plaq = Plaq + sitePlaq; + } + } + } + + // sum_mu_nu_x Re tr U_mu_nu(x) + static RealD sumPlaquette(const GaugeLorentz &Umu) { + std::vector U(Nd, Umu.Grid()); + for (int mu = 0; mu < Nd; mu++) { + // Umu is a GaugeLorentz object, and as such has a non-trivial Lorentz index. We can + // access the element in the mu Lorentz index with this PeekIndex syntax. + U[mu] = PeekIndex(Umu, mu); + } + ComplexField Plaq(Umu.Grid()); + sitePlaquette(Plaq, U); + // I guess this should be the line that sums over all space-time sites. + auto Tp = sum(Plaq); + // Until now, we have been working with objects inside the tensor nest. This TensorRemove gets + // rid of the tensor nest to return whatever is inside. + auto p = TensorRemove(Tp); + return p.real(); + } + + // < Re tr U_mu_nu(x) > + static RealD avgPlaquette(const GaugeLorentz &Umu) { + // Real double type + RealD sumplaq = sumPlaquette(Umu); + // gSites() is the number of global sites. there is also lSites() for local sites. + double vol = Umu.Grid()->gSites(); + // The number of orientations. 4*3/2=6 for Nd=4, as known. + double faces = (1.0 * Nd * (Nd - 1)) / 2.0; + return sumplaq / vol / faces / Nc; + } +}; + + +// Next we show an example of how to construct an input parameter class. We first inherit +// from Serializable. Then all class data members have to be defined using the +// GRID_SERIALIZABLE_CLASS_MEMBERS macro. This variadic macro allows for arbitrarily many +// class data members. In the below case, we make a parameter file holding the configuration +// name. Here, it expects the name to be labeled with "conf_name" in the configuration file. +struct ConfParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS( + ConfParameters, + std::string, conf_name); + + template + ConfParameters(Reader& Reader){ + // If we are reading an XML file, it should be structured like: + // + // + // l20t20b06498a_nersc.302500 + // + // + read(Reader, "parameters", *this); + } +}; + + + +// This syntax lets you pass command line arguments to main. An asterisk means that what follows is +// a pointer. Two asterisks means what follows is a pointer to an array. +int main (int argc, char **argv) +{ + // This initializes Grid. Some command line options include + // --mpi n.n.n.n + // --threads n + // --grid n.n.n.n + Grid_init(&argc, &argv); + + // This is where you would specify a custom lattice size, if not from the command line. Here + // Nd is a global quantity that is currently set to 4. + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + Coordinate latt_size = GridDefaultLatt(); + + // Instantiate the spacetime Grid on which everything will be built. + GridCartesian GRID(latt_size,simd_layout,mpi_layout); + + // The PeriodicGimplD type is what you want for gauge matrices. There is also a LatticeGaugeFieldD + // type that you can use, which will work perfectly with what follows. + PeriodicGimplD::Field U(&GRID); + + // Here we read in the parameter file params.json to get conf_name. The last argument is what the + // top organizational level is called in the param file. + XmlReader Reader("Example_plaquette.xml",false, "grid"); + ConfParameters param(Reader); + + // Load a lattice from SIMULATeQCD into U. SIMULATeQCD finds plaquette = 0.6381995717 + FieldMetaData header; + NerscIO::readConfiguration(U, header, param.conf_name); + + // Let's see what we find. + RealD plaq = WLoops::avgPlaquette(U); + + // This is how you make log messages. + std::cout << GridLogMessage << std::setprecision(std::numeric_limits::digits10 + 1) << "Plaquette = " << plaq << std::endl; + + // To wrap things up. + Grid_finalize(); +} \ No newline at end of file diff --git a/tests/debug/Test_padded_cell.cc b/tests/debug/Test_padded_cell.cc index 4fb461fe..f110df46 100644 --- a/tests/debug/Test_padded_cell.cc +++ b/tests/debug/Test_padded_cell.cc @@ -32,6 +32,7 @@ Author: Peter Boyle using namespace std; using namespace Grid; +// This is to optimize the SIMD template void gpermute(vobj & inout,int perm){ vobj tmp=inout; if (perm & 0x1 ) { permute(inout,tmp,0); tmp=inout;} @@ -39,7 +40,8 @@ template void gpermute(vobj & inout,int perm){ if (perm & 0x4 ) { permute(inout,tmp,2); tmp=inout;} if (perm & 0x8 ) { permute(inout,tmp,3); tmp=inout;} } - + + int main (int argc, char ** argv) { Grid_init(&argc,&argv); @@ -47,20 +49,21 @@ int main (int argc, char ** argv) Coordinate latt_size = GridDefaultLatt(); Coordinate simd_layout= GridDefaultSimd(Nd,vComplexD::Nsimd()); Coordinate mpi_layout = GridDefaultMpi(); - std::cout << " mpi "<({45,12,81,9})); LatticeGaugeField Umu(&GRID); - + pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); SU::HotConfiguration(pRNG,Umu); Real plaq=WilsonLoops::avgPlaquette(Umu); LatticeComplex trplaq(&GRID); + // Store Umu in U. Peek/Poke mean respectively getElement/setElement. std::vector U(Nd, Umu.Grid()); for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(Umu, mu); @@ -70,9 +73,7 @@ int main (int argc, char ** argv) LatticeComplex cplaq(&GRID); cplaq=Zero(); - ///////////////////////////////////////////////// // Create a padded cell of extra padding depth=1 - ///////////////////////////////////////////////// int depth = 1; PaddedCell Ghost(depth,&GRID); LatticeGaugeField Ughost = Ghost.Exchange(Umu); @@ -114,18 +115,25 @@ int main (int argc, char ** argv) } #endif - ///// Array for the site plaquette + // Array for the site plaquette GridBase *GhostGrid = Ughost.Grid(); LatticeComplex gplaq(GhostGrid); - + + // Now we're going to put together the "stencil" that will be useful to us when + // calculating the plaquette. Our eventual goal is to make the product + // Umu(x) Unu(x+mu) Umu^dag(x+nu) Unu^dag(x), + // which requires, in order, the sites x, x+mu, x+nu, and x. We arrive at these + // sites relative to x through "shifts", which is represented here by a 4-d + // vector of 0s (no movement) and 1s (shift one unit) at each site. The + // "stencil" is the set of all these shifts. std::vector shifts; for(int mu=0;mu_offset; - int o1 = SE1->_offset; - int o2 = SE2->_offset; - int o3 = SE3->_offset; - - auto U0 = U_v[o0](mu); - auto U1 = U_v[o1](nu); - auto U2 = adj(U_v[o2](mu)); - auto U3 = adj(U_v[o3](nu)); + // Before doing accelerator stuff, there is an opening and closing of "Views". I guess the + // "Views" are stored in *_v variables listed below. + autoView( gp_v , gplaq, CpuWrite); + autoView( t_v , trplaq, CpuRead); + autoView( U_v , Ughost, CpuRead); - gpermute(U0,SE0->_permute); - gpermute(U1,SE1->_permute); - gpermute(U2,SE2->_permute); - gpermute(U3,SE3->_permute); - - gp_v[ss]() =gp_v[ss]() + trace( U0*U1*U2*U3 ); - s=s+4; - } - } + // This is now a loop over stencil shift elements. That is, s increases as we make our + // way through the spacetimes sites, but also as we make our way around the plaquette. + for(int ss=0;ss_offset; + int o1 = SE1->_offset; + int o2 = SE2->_offset; + int o3 = SE3->_offset; + + auto U0 = U_v[o0](mu); + auto U1 = U_v[o1](nu); + auto U2 = adj(U_v[o2](mu)); + auto U3 = adj(U_v[o3](nu)); + + gpermute(U0,SE0->_permute); + gpermute(U1,SE1->_permute); + gpermute(U2,SE2->_permute); + gpermute(U3,SE3->_permute); + + gp_v[ss]() =gp_v[ss]() + trace( U0*U1*U2*U3 ); + s=s+4; + } } } + + // Here is my understanding of this part: The padded cell has its own periodic BCs, so + // if I take a step to the right at the right-most side of the cell, I end up on the + // left-most side. This means that the plaquettes in the padding are wrong. Luckily + // all we care about are the plaquettes in the cell, which we obtain from Extract. cplaq = Ghost.Extract(gplaq); RealD vol = cplaq.Grid()->gSites(); RealD faces = (Nd * (Nd-1))/2; diff --git a/tests/smearing/Test_fatLinks.cc b/tests/smearing/Test_fatLinks.cc new file mode 100644 index 00000000..04d5b165 --- /dev/null +++ b/tests/smearing/Test_fatLinks.cc @@ -0,0 +1,181 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/smearing/Test_fatLinks.cc + +Copyright (C) 2023 + +Author: D. A. Clarke + +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 +*************************************************************************************/ +/* + @file Test_fatLinks.cc + @brief test of the HISQ smearing +*/ + + +#include +#include +#include +#include +using namespace Grid; + + +/*! @brief parameter file to easily adjust Nloop */ +struct ConfParameters: Serializable { + GRID_SERIALIZABLE_CLASS_MEMBERS( + ConfParameters, + int, benchmark, + int, Nloop); + + template + ConfParameters(Reader& Reader){ + read(Reader, "parameters", *this); + } +}; + + +bool testSmear(GridCartesian& GRID, LatticeGaugeFieldD Umu, LatticeGaugeFieldD Usmr, LatticeGaugeFieldD Unaik, + LatticeGaugeFieldD Ucontrol, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) { + Smear_HISQ hisq_fat(&GRID,c1,cnaik,c3,c5,c7,clp); + LatticeGaugeFieldD diff(&GRID), Uproj(&GRID); + hisq_fat.smear(Usmr, Unaik, Umu); + bool result; + if (cnaik < 1e-30) { // Testing anything but Naik term + diff = Ucontrol-Usmr; + auto absDiff = norm2(diff)/norm2(Ucontrol); + if (absDiff < 1e-30) { + Grid_pass(" |Umu-Usmr|/|Umu| = ",absDiff); + result = true; + } else { + Grid_error(" |Umu-Usmr|/|Umu| = ",absDiff); + result = false; + } + } else { // Testing Naik specifically + diff = Ucontrol-Unaik; + auto absDiff = norm2(diff)/norm2(Ucontrol); + if (absDiff < 1e-30) { + Grid_pass(" |Umu-Unaik|/|Umu| = ",absDiff); + result = true; + } else { + Grid_error(" |Umu-Unaik|/|Umu| = ",absDiff); + result = false; + } + hisq_fat.projectU3(Uproj,Ucontrol); +// NerscIO::writeConfiguration(Unaik,"nersc.l8t4b3360.naik"); + } + return result; +} + + +int main (int argc, char** argv) { + + // Params for the test. + int Ns = 8; + int Nt = 4; + Coordinate latt_size(Nd,0); latt_size[0]=Ns; latt_size[1]=Ns; latt_size[2]=Ns; latt_size[3]=Nt; + std::string conf_in = "nersc.l8t4b3360"; + int threads = GridThread::GetThreads(); + + typedef LatticeGaugeFieldD LGF; + + // Initialize the Grid + Grid_init(&argc,&argv); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + Grid_log("mpi = ",mpi_layout); + Grid_log("simd = ",simd_layout); + Grid_log("latt = ",latt_size); + Grid_log("threads = ",threads); + GridCartesian GRID(latt_size,simd_layout,mpi_layout); + + XmlReader Reader("fatParams.xml",false,"grid"); + ConfParameters param(Reader); + if(param.benchmark) Grid_log(" Nloop = ",param.Nloop); + + LGF Umu(&GRID), Usmr(&GRID), Unaik(&GRID), Ucontrol(&GRID); + + // Read the configuration into Umu + FieldMetaData header; + NerscIO::readConfiguration(Umu, header, conf_in); + + bool pass=true; + + // Carry out various tests + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357lplink.control"); + pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,-1/8.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.357link.control"); + pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,1/384.,0.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.35link.control"); + pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,1/64.,0.,0.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.3link.control"); + pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,1/8.,0.,1/16.,0.,0.,0.); + NerscIO::readConfiguration(Ucontrol, header, "nersc.l8t4b3360.naik.control"); + pass *= testSmear(GRID,Umu,Usmr,Unaik,Ucontrol,0.,0.8675309,0.,0.,0.,0.); + + if(pass){ + Grid_pass("All tests passed."); + } else { + Grid_error("At least one test failed."); + } + + // Test a C-style instantiation + double path_coeff[6] = {1, 2, 3, 4, 5, 6}; + Smear_HISQ hisq_fat_Cstyle(&GRID,path_coeff); + + if (param.benchmark) { + + autoView(U_v, Umu, CpuRead); // Gauge accessor + + // Read in lattice sequentially, Nloop times + double lookupTime = 0.; + for(int i=0;i