From 831ca4e3bf8e0b4231f395b2af9308e007b73186 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Tue, 14 Mar 2017 14:55:18 +0900 Subject: [PATCH 1/8] Added Scalar action for fields in the adjoint representation --- lib/qcd/action/Actions.h | 5 + lib/qcd/action/scalar/ScalarAction.h | 61 ++++++----- lib/qcd/action/scalar/ScalarImpl.h | 93 ++++++++-------- .../action/scalar/ScalarInteractionAction.h | 84 +++++++-------- lib/qcd/hmc/GenericHMCrunner.h | 3 + lib/qcd/representations/hmc_types.h | 2 +- tests/hmc/Test_hmc_ScalarActionNxN.cc | 100 ++++++++++++++++++ 7 files changed, 227 insertions(+), 121 deletions(-) create mode 100644 tests/hmc/Test_hmc_ScalarActionNxN.cc diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index daf64f3d..0214b8f4 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -69,6 +69,7 @@ Author: paboyle //////////////////////////////////////////// #include #include +#include namespace Grid { namespace QCD { @@ -106,6 +107,10 @@ typedef ScalarAction ScalarActionR; typedef ScalarAction ScalarActionF; typedef ScalarAction ScalarActionD; +typedef ScalarInteractionAction ScalarAdjActionR; +typedef ScalarInteractionAction ScalarAdjActionF; +typedef ScalarInteractionAction ScalarAdjActionD; + }} //////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/lib/qcd/action/scalar/ScalarAction.h b/lib/qcd/action/scalar/ScalarAction.h index f10ec9a6..2c82d2e3 100644 --- a/lib/qcd/action/scalar/ScalarAction.h +++ b/lib/qcd/action/scalar/ScalarAction.h @@ -6,10 +6,10 @@ Copyright (C) 2015 -Author: Azusa Yamaguchi -Author: Peter Boyle -Author: neo -Author: paboyle + Author: Azusa Yamaguchi + Author: Peter Boyle + Author: neo + Author: paboyle 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 @@ -35,50 +35,49 @@ directory namespace Grid { // FIXME drop the QCD namespace everywhere here - - template - class ScalarAction : public QCD::Action { - public: + +template +class ScalarAction : public QCD::Action { + public: INHERIT_FIELD_TYPES(Impl); - - private: + + private: RealD mass_square; RealD lambda; - - public: - ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){}; - virtual std::string LogParameters(){ + public: + ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l) {} + + virtual std::string LogParameters() { std::stringstream sstream; sstream << GridLogMessage << "[ScalarAction] lambda : " << lambda << std::endl; sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl; return sstream.str(); - } - - virtual std::string action_name(){return "ScalarAction";} - - virtual void refresh(const Field &U, - GridParallelRNG &pRNG){}; // noop as no pseudoferms - + virtual std::string action_name() {return "ScalarAction";} + + virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} // noop as no pseudoferms + virtual RealD S(const Field &p) { return (mass_square * 0.5 + QCD::Nd) * ScalarObs::sumphisquared(p) + - (lambda / 24.) * ScalarObs::sumphifourth(p) + - ScalarObs::sumphider(p); + (lambda / 24.) * ScalarObs::sumphifourth(p) + + ScalarObs::sumphider(p); }; - + virtual void deriv(const Field &p, - Field &force) { + Field &force) { Field tmp(p._grid); Field p2(p._grid); ScalarObs::phisquared(p2, p); tmp = -(Cshift(p, 0, -1) + Cshift(p, 0, 1)); for (int mu = 1; mu < QCD::Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1); - - force=+(mass_square + 2. * QCD::Nd) * p + (lambda / 6.) * p2 * p + tmp; - }; - }; - -} // Grid + + force =+(mass_square + 2. * QCD::Nd) * p + (lambda / 6.) * p2 * p + tmp; + } +}; + + + +} // namespace Grid #endif // SCALAR_ACTION_H diff --git a/lib/qcd/action/scalar/ScalarImpl.h b/lib/qcd/action/scalar/ScalarImpl.h index ee2d2fb8..6d14b61a 100644 --- a/lib/qcd/action/scalar/ScalarImpl.h +++ b/lib/qcd/action/scalar/ScalarImpl.h @@ -5,96 +5,99 @@ namespace Grid { //namespace QCD { - template - class ScalarImplTypes { - public: +template +class ScalarImplTypes { + public: typedef S Simd; - + template using iImplField = iScalar > >; - + typedef iImplField SiteField; - - + typedef Lattice Field; - - static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){ + + static inline void generate_momenta(Field& P, GridParallelRNG& pRNG) { gaussian(pRNG, P); } - + static inline Field projectForce(Field& P){return P;} - - static inline void update_field(Field& P, Field& U, double ep){ + + static inline void update_field(Field& P, Field& U, double ep) { U += P*ep; } - - static inline RealD FieldSquareNorm(Field& U){ + + static inline RealD FieldSquareNorm(Field& U) { return (- sum(trace(U*U))/2.0); } - + static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) { gaussian(pRNG, U); } - + static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) { gaussian(pRNG, U); } - + static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) { U = 1.0; } - + }; template - class ScalarMatrixImplTypes { + class ScalarAdjMatrixImplTypes { public: typedef S Simd; - template using iImplField = iScalar > >; - + typedef iImplField SiteField; - - + typedef Lattice Field; - - static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){ - gaussian(pRNG, P); + + static inline void generate_momenta(Field& P, GridParallelRNG& pRNG) { + QCD::SU::GaussianFundamentalLieAlgebraMatrix(pRNG, P); } - - static inline Field projectForce(Field& P){return P;} - - static inline void update_field(Field& P, Field& U, double ep){ + + static inline Field projectForce(Field& P) {return P;} + + static inline void update_field(Field& P, Field& U, double ep) { U += P*ep; } - - static inline RealD FieldSquareNorm(Field& U){ - return (TensorRemove(- sum(trace(U*U))*0.5).real()); + + static inline RealD FieldSquareNorm(Field& U) { + return (TensorRemove(sum(trace(U*U))).real()); } - + static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) { - gaussian(pRNG, U); + QCD::SU::LieRandomize(pRNG, U); } - + static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) { - gaussian(pRNG, U); + QCD::SU::LieRandomize(pRNG, U, 0.01); } - + static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) { - U = 1.0; + U = zero; } - + }; - - + + typedef ScalarImplTypes ScalarImplR; typedef ScalarImplTypes ScalarImplF; typedef ScalarImplTypes ScalarImplD; - - //} -} + + // Hardcoding here the size of the matrices + typedef ScalarAdjMatrixImplTypes ScalarAdjImplR; + typedef ScalarAdjMatrixImplTypes ScalarAdjImplF; + typedef ScalarAdjMatrixImplTypes ScalarAdjImplD; + + + //} +} #endif diff --git a/lib/qcd/action/scalar/ScalarInteractionAction.h b/lib/qcd/action/scalar/ScalarInteractionAction.h index bd54a010..2607b041 100644 --- a/lib/qcd/action/scalar/ScalarInteractionAction.h +++ b/lib/qcd/action/scalar/ScalarInteractionAction.h @@ -6,10 +6,7 @@ Copyright (C) 2015 -Author: Azusa Yamaguchi -Author: Peter Boyle -Author: neo -Author: paboyle + Author: Guido Cossu 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 @@ -30,55 +27,54 @@ directory *************************************************************************************/ /* END LEGAL */ -#ifndef SCALAR_ACTION_H -#define SCALAR_ACTION_H +#ifndef SCALAR_INT_ACTION_H +#define SCALAR_INT_ACTION_H namespace Grid { // FIXME drop the QCD namespace everywhere here - - template - class ScalarInteractionAction : public QCD::Action { - public: - INHERIT_FIELD_TYPES(Impl); - - private: + +template +class ScalarInteractionAction : public QCD::Action { RealD mass_square; RealD lambda; - - public: - ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){}; - virtual std::string LogParameters(){ + public: + INHERIT_FIELD_TYPES(Impl); + ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l) {} + + virtual std::string LogParameters() { std::stringstream sstream; sstream << GridLogMessage << "[ScalarAction] lambda : " << lambda << std::endl; sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl; return sstream.str(); - } - - virtual std::string action_name(){return "ScalarAction";} - - virtual void refresh(const Field &U, - GridParallelRNG &pRNG){}; // noop as no pseudoferms - - virtual RealD S(const Field &p) { - return (mass_square * 0.5 + QCD::Nd) * ScalarObs::sumphisquared(p) + - (lambda / 24.) * ScalarObs::sumphifourth(p) + - ScalarObs::sumphider(p); - }; - - virtual void deriv(const Field &p, - Field &force) { - Field tmp(p._grid); - Field p2(p._grid); - ScalarObs::phisquared(p2, p); - tmp = -(Cshift(p, 0, -1) + Cshift(p, 0, 1)); - for (int mu = 1; mu < QCD::Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1); - - force=+(mass_square + 2. * QCD::Nd) * p + (lambda / 6.) * p2 * p + tmp; - }; - }; - -} // Grid -#endif // SCALAR_ACTION_H + virtual std::string action_name() {return "ScalarAction";} + + virtual void refresh(const Field &U, + GridParallelRNG &pRNG) {} // noop as no pseudoferms + + virtual RealD S(const Field &p) { + Field action(p._grid); + Field pshift(p._grid); + Field phisquared(p._grid); + phisquared = p*p; + action = (2.0*QCD::Nd + mass_square)*phisquared + lambda*phisquared*phisquared; + for (int mu = 0; mu < QCD::Nd; mu++) { + pshift = Cshift(p, mu, +1); // not efficient implement with stencils + action -= pshift*p + p*pshift; + } + return -(TensorRemove(sum(trace(action)))).real(); + }; + + virtual void deriv(const Field &p, + Field &force) { + force = (2.0*QCD::Nd + mass_square)*p + 2.0*lambda*p*p*p; + // following is inefficient + for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); + } +}; + +} // namespace Grid + +#endif // SCALAR_INT_ACTION_H diff --git a/lib/qcd/hmc/GenericHMCrunner.h b/lib/qcd/hmc/GenericHMCrunner.h index 66b16435..a97fb4e4 100644 --- a/lib/qcd/hmc/GenericHMCrunner.h +++ b/lib/qcd/hmc/GenericHMCrunner.h @@ -202,6 +202,9 @@ using GenericHMCRunnerTemplate = HMCWrapperTemplate ScalarGenericHMCRunner; +typedef HMCWrapperTemplate + ScalarAdjGenericHMCRunner; + } // namespace QCD } // namespace Grid diff --git a/lib/qcd/representations/hmc_types.h b/lib/qcd/representations/hmc_types.h index 3701c9b2..b4991941 100644 --- a/lib/qcd/representations/hmc_types.h +++ b/lib/qcd/representations/hmc_types.h @@ -62,7 +62,7 @@ class Representations { typedef Representations NoHirep; typedef Representations > ScalarFields; - //typedef Representations > ScalarMatrixFields; +typedef Representations > ScalarMatrixFields; // Helper classes to access the elements // Strips the first N parameters from the tuple diff --git a/tests/hmc/Test_hmc_ScalarActionNxN.cc b/tests/hmc/Test_hmc_ScalarActionNxN.cc new file mode 100644 index 00000000..8b93efde --- /dev/null +++ b/tests/hmc/Test_hmc_ScalarActionNxN.cc @@ -0,0 +1,100 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_hmc_WilsonFermionGauge.cc + +Copyright (C) 2016 + +Author: Guido Cossu + +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 +*************************************************************************************/ +/* END LEGAL */ +#include +namespace Grid{ +class ScalarActionParameters : Serializable { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarActionParameters, + double, mass_squared, + double, lambda); +}; + +} +int main(int argc, char **argv) { + using namespace Grid; + using namespace Grid::QCD; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef ScalarAdjGenericHMCRunner HMCWrapper; // Uses the default minimum norm, real scalar fields + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + HMCWrapper TheHMC; + + // Grid from the command line + GridModule ScalarGrid; + ScalarGrid.set_full( SpaceTimeGrid::makeFourDimGrid( + GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), + GridDefaultMpi())); + ScalarGrid.set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(ScalarGrid.get_full())); + TheHMC.Resources.AddGrid("scalar", ScalarGrid); + // Possibile to create the module by hand + // hardcoding parameters or using a Reader + + // Checkpointer definition + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_scalar_lat"; + CPparams.rng_prefix = "ckpoint_scalar_rng"; + CPparams.saveInterval = 50; + CPparams.format = "IEEE64BIG"; + + TheHMC.Resources.LoadBinaryCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + ///////////////////////////////////////////////////////////// + // Collect actions, here use more encapsulation + + // Scalar action in adjoint representation + ScalarActionParameters SPar; + SPar.mass_squared = 0.5; + SPar.lambda = 0.1; + ScalarAdjActionR Saction(SPar.mass_squared, SPar.lambda); + + // Collect actions + ActionLevel Level1(1); + Level1.push_back(&Saction); + TheHMC.TheAction.push_back(Level1); + ///////////////////////////////////////////////////////////// + + // HMC parameters are serialisable + TheHMC.Parameters.MD.MDsteps = 10; + TheHMC.Parameters.MD.trajL = 1.0; + + TheHMC.ReadCommandLine(argc, argv); + TheHMC.Run(); + + Grid_finalize(); + +} // main From 38806343a873ea10264c79103db31182d6770947 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Wed, 15 Mar 2017 15:16:16 +0900 Subject: [PATCH 2/8] Improving efficiency of the force term --- .../action/scalar/ScalarInteractionAction.h | 91 ++++++++++++++++--- tests/Test_stencil.cc | 43 +++++---- tests/hmc/Test_hmc_ScalarActionNxN.cc | 11 +-- 3 files changed, 104 insertions(+), 41 deletions(-) diff --git a/lib/qcd/action/scalar/ScalarInteractionAction.h b/lib/qcd/action/scalar/ScalarInteractionAction.h index 2607b041..5a322a5e 100644 --- a/lib/qcd/action/scalar/ScalarInteractionAction.h +++ b/lib/qcd/action/scalar/ScalarInteractionAction.h @@ -30,17 +30,34 @@ directory #ifndef SCALAR_INT_ACTION_H #define SCALAR_INT_ACTION_H + +// Note: this action can completely absorb the ScalarAction for real float fields +// use the scalarObjs to generalise the structure + namespace Grid { // FIXME drop the QCD namespace everywhere here template class ScalarInteractionAction : public QCD::Action { +public: + INHERIT_FIELD_TYPES(Impl); +private: RealD mass_square; RealD lambda; + + typedef typename Field::vector_object vobj; + typedef CartesianStencil Stencil; + + SimpleCompressor compressor; + int npoint = 8; + std::vector directions = {0,1,2,3,0,1,2,3}; // forcing 4 dimensions + std::vector displacements = {1,1,1,1, -1,-1,-1,-1}; + + public: - INHERIT_FIELD_TYPES(Impl); - ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l) {} + + ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l){} virtual std::string LogParameters() { std::stringstream sstream; @@ -51,27 +68,75 @@ class ScalarInteractionAction : public QCD::Action { virtual std::string action_name() {return "ScalarAction";} - virtual void refresh(const Field &U, - GridParallelRNG &pRNG) {} // noop as no pseudoferms + virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} virtual RealD S(const Field &p) { - Field action(p._grid); - Field pshift(p._grid); - Field phisquared(p._grid); + static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); + phiStencil.HaloExchange(p, compressor); + + Field action(p._grid), pshift(p._grid), phisquared(p._grid); phisquared = p*p; action = (2.0*QCD::Nd + mass_square)*phisquared + lambda*phisquared*phisquared; for (int mu = 0; mu < QCD::Nd; mu++) { - pshift = Cshift(p, mu, +1); // not efficient implement with stencils - action -= pshift*p + p*pshift; + // pshift = Cshift(p, mu, +1); // not efficient, implement with stencils + PARALLEL_FOR_LOOP + for (int i = 0; i < p._grid->oSites(); i++) { + int permute_type; + StencilEntry *SE; + vobj temp2; + vobj *temp; + vobj *t_p; + + SE = phiStencil.GetEntry(permute_type, mu, i); + t_p = &p._odata[i]; + if ( SE->_is_local ) { + temp = &p._odata[SE->_offset]; + if ( SE->_permute ) { + permute(temp2, *temp, permute_type); + action._odata[i] -= temp2*(*t_p) + (*t_p)*temp2; + } else { + action._odata[i] -= *temp*(*t_p) + (*t_p)*(*temp); + } + } else { + action._odata[i] -= phiStencil.CommBuf()[SE->_offset]*(*t_p) + (*t_p)*phiStencil.CommBuf()[SE->_offset]; + } + } + // action -= pshift*p + p*pshift; } + // NB the trace in the algebra is normalised to 1/2 + // minus sign coming from the antihermitian fields return -(TensorRemove(sum(trace(action)))).real(); }; - virtual void deriv(const Field &p, - Field &force) { + virtual void deriv(const Field &p, Field &force) { force = (2.0*QCD::Nd + mass_square)*p + 2.0*lambda*p*p*p; - // following is inefficient - for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); + // move this outside + static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); + phiStencil.HaloExchange(p, compressor); + + //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); + for (int point = 0; point < npoint; point++) { + PARALLEL_FOR_LOOP + for (int i = 0; i < p._grid->oSites(); i++) { + vobj *temp; + vobj temp2; + int permute_type; + StencilEntry *SE; + SE = phiStencil.GetEntry(permute_type, point, i); + + if ( SE->_is_local ) { + temp = &p._odata[SE->_offset]; + if ( SE->_permute ) { + permute(temp2, *temp, permute_type); + force._odata[i] -= temp2; + } else { + force._odata[i] -= *temp; + } + } else { + force._odata[i] -= phiStencil.CommBuf()[SE->_offset]; + } + } + } } }; diff --git a/tests/Test_stencil.cc b/tests/Test_stencil.cc index 1b71b8a5..1d35e1bb 100644 --- a/tests/Test_stencil.cc +++ b/tests/Test_stencil.cc @@ -1,6 +1,6 @@ /************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid + Grid physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_stencil.cc @@ -33,9 +33,8 @@ using namespace std; using namespace Grid; using namespace Grid::QCD; -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); +int main(int argc, char ** argv) { + Grid_init(&argc, &argv); // typedef LatticeColourMatrix Field; typedef LatticeComplex Field; @@ -47,7 +46,7 @@ int main (int argc, char ** argv) std::vector mpi_layout = GridDefaultMpi(); double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3]; - + GridCartesian Fine(latt_size,simd_layout,mpi_layout); GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout); GridParallelRNG fRNG(&Fine); @@ -55,14 +54,14 @@ int main (int argc, char ** argv) // fRNG.SeedRandomDevice(); std::vector seeds({1,2,3,4}); fRNG.SeedFixedIntegers(seeds); - + Field Foo(&Fine); Field Bar(&Fine); Field Check(&Fine); Field Diff(&Fine); LatticeComplex lex(&Fine); - lex = zero; + lex = zero; random(fRNG,Foo); gaussian(fRNG,Bar); @@ -98,7 +97,7 @@ int main (int argc, char ** argv) Fine.oCoorFromOindex(ocoor,o); ocoor[dir]=(ocoor[dir]+disp)%Fine._rdimensions[dir]; } - + SimpleCompressor compress; myStencil.HaloExchange(Foo,compress); @@ -106,16 +105,16 @@ int main (int argc, char ** argv) // Implement a stencil code that should agree with cshift! for(int i=0;ioSites();i++){ - + int permute_type; StencilEntry *SE; SE = myStencil.GetEntry(permute_type,0,i); - + if ( SE->_is_local && SE->_permute ) permute(Check._odata[i],Foo._odata[SE->_offset],permute_type); else if (SE->_is_local) Check._odata[i] = Foo._odata[SE->_offset]; - else + else Check._odata[i] = myStencil.CommBuf()[SE->_offset]; } @@ -144,7 +143,7 @@ int main (int argc, char ** argv) <<") " < compress; EStencil.HaloExchange(EFoo,compress); OStencil.HaloExchange(OFoo,compress); - + Bar = Cshift(Foo,dir,disp); if ( disp & 0x1 ) { ECheck.checkerboard = Even; OCheck.checkerboard = Odd; - } else { + } else { ECheck.checkerboard = Odd; OCheck.checkerboard = Even; } @@ -206,7 +205,7 @@ int main (int argc, char ** argv) permute(OCheck._odata[i],EFoo._odata[SE->_offset],permute_type); else if (SE->_is_local) OCheck._odata[i] = EFoo._odata[SE->_offset]; - else + else OCheck._odata[i] = EStencil.CommBuf()[SE->_offset]; } for(int i=0;ioSites();i++){ @@ -214,18 +213,18 @@ int main (int argc, char ** argv) StencilEntry *SE; SE = OStencil.GetEntry(permute_type,0,i); // std::cout << "ODD source "<< i<<" -> " <_offset << " "<< SE->_is_local<_is_local && SE->_permute ) permute(ECheck._odata[i],OFoo._odata[SE->_offset],permute_type); else if (SE->_is_local) ECheck._odata[i] = OFoo._odata[SE->_offset]; - else + else ECheck._odata[i] = OStencil.CommBuf()[SE->_offset]; } - + setCheckerboard(Check,ECheck); setCheckerboard(Check,OCheck); - + Real nrmC = norm2(Check); Real nrmB = norm2(Bar); Diff = Check-Bar; @@ -248,10 +247,10 @@ int main (int argc, char ** argv) diff =norm2(ddiff); if ( diff > 0){ std::cout <<"Coor (" << coor[0]<<","< -namespace Grid{ +namespace Grid { class ScalarActionParameters : Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarActionParameters, @@ -44,7 +44,7 @@ int main(int argc, char **argv) { // here make a routine to print all the relevant information on the run std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; - // Typedefs to simplify notation + // Typedefs to simplify notation typedef ScalarAdjGenericHMCRunner HMCWrapper; // Uses the default minimum norm, real scalar fields //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -52,7 +52,7 @@ int main(int argc, char **argv) { // Grid from the command line GridModule ScalarGrid; - ScalarGrid.set_full( SpaceTimeGrid::makeFourDimGrid( + ScalarGrid.set_full(SpaceTimeGrid::makeFourDimGrid( GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi())); ScalarGrid.set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(ScalarGrid.get_full())); @@ -89,12 +89,11 @@ int main(int argc, char **argv) { ///////////////////////////////////////////////////////////// // HMC parameters are serialisable - TheHMC.Parameters.MD.MDsteps = 10; + TheHMC.Parameters.MD.MDsteps = 20; TheHMC.Parameters.MD.trajL = 1.0; TheHMC.ReadCommandLine(argc, argv); TheHMC.Run(); Grid_finalize(); - -} // main +} // main From 038b6ee9cdfc5902b27a8645b1f1758c9db3656f Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 16 Mar 2017 01:09:24 +0900 Subject: [PATCH 3/8] Fixing JSON compilation error --- lib/json/json.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/json/json.hpp b/lib/json/json.hpp index 97214f0b..bfb38c3e 100644 --- a/lib/json/json.hpp +++ b/lib/json/json.hpp @@ -64,7 +64,7 @@ SOFTWARE. #endif #elif defined(__GNUC__) #define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) - #if GCC_VERSION < 40900 + #if GCC_VERSION < 40800 #error "unsupported GCC version - see https://github.com/nlohmann/json#supported-compilers" #endif #endif From 7b03d8d0879d7f7922b8867eefa9346cb0e5c425 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Wed, 5 Apr 2017 16:17:46 +0100 Subject: [PATCH 4/8] Fixing the remaining merge conflicts --- lib/qcd/action/scalar/Scalar.h | 5 +++++ tests/Test_stencil.cc | 7 ------- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/lib/qcd/action/scalar/Scalar.h b/lib/qcd/action/scalar/Scalar.h index e5bea275..cae38360 100644 --- a/lib/qcd/action/scalar/Scalar.h +++ b/lib/qcd/action/scalar/Scalar.h @@ -31,6 +31,7 @@ directory #include #include +#include namespace Grid { namespace QCD { @@ -39,6 +40,10 @@ namespace QCD { typedef ScalarAction ScalarActionF; typedef ScalarAction ScalarActionD; + typedef ScalarInteractionAction ScalarAdjActionR; + typedef ScalarInteractionAction ScalarAdjActionF; + typedef ScalarInteractionAction ScalarAdjActionD; + } } diff --git a/tests/Test_stencil.cc b/tests/Test_stencil.cc index 2a4744f3..fa4b0b57 100644 --- a/tests/Test_stencil.cc +++ b/tests/Test_stencil.cc @@ -189,13 +189,6 @@ int main(int argc, char ** argv) { SimpleCompressor compress; -<<<<<<< HEAD - EStencil.HaloExchange(EFoo,compress); - OStencil.HaloExchange(OFoo,compress); - -======= - ->>>>>>> feature/hmc_generalise Bar = Cshift(Foo,dir,disp); if ( disp & 0x1 ) { From 741bc836f69d37623cba76cf4aee06dee3f6c84e Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 5 May 2017 17:36:43 +0100 Subject: [PATCH 5/8] Exposing support for Ncolours and Ndimensions and JSON input file for the ScalarAction --- lib/qcd/action/scalar/Scalar.h | 6 +- lib/qcd/action/scalar/ScalarImpl.h | 5 +- .../action/scalar/ScalarInteractionAction.h | 152 +++++++++--------- lib/qcd/hmc/GenericHMCrunner.h | 3 + lib/qcd/representations/hmc_types.h | 3 + lib/stencil/Stencil.h | 2 +- tests/hmc/Test_hmc_ScalarActionNxN.cc | 104 ++++++++---- 7 files changed, 168 insertions(+), 107 deletions(-) diff --git a/lib/qcd/action/scalar/Scalar.h b/lib/qcd/action/scalar/Scalar.h index cae38360..485a6765 100644 --- a/lib/qcd/action/scalar/Scalar.h +++ b/lib/qcd/action/scalar/Scalar.h @@ -40,9 +40,9 @@ namespace QCD { typedef ScalarAction ScalarActionF; typedef ScalarAction ScalarActionD; - typedef ScalarInteractionAction ScalarAdjActionR; - typedef ScalarInteractionAction ScalarAdjActionF; - typedef ScalarInteractionAction ScalarAdjActionD; + template using ScalarAdjActionR = ScalarInteractionAction, Dimensions>; + template using ScalarAdjActionF = ScalarInteractionAction, Dimensions>; + template using ScalarAdjActionD = ScalarInteractionAction, Dimensions>; } } diff --git a/lib/qcd/action/scalar/ScalarImpl.h b/lib/qcd/action/scalar/ScalarImpl.h index 6d14b61a..8b5e3aa2 100644 --- a/lib/qcd/action/scalar/ScalarImpl.h +++ b/lib/qcd/action/scalar/ScalarImpl.h @@ -96,7 +96,10 @@ class ScalarImplTypes { typedef ScalarAdjMatrixImplTypes ScalarAdjImplF; typedef ScalarAdjMatrixImplTypes ScalarAdjImplD; - + template using ScalarNxNAdjImplR = ScalarAdjMatrixImplTypes; + template using ScalarNxNAdjImplF = ScalarAdjMatrixImplTypes; + template using ScalarNxNAdjImplD = ScalarAdjMatrixImplTypes; + //} } diff --git a/lib/qcd/action/scalar/ScalarInteractionAction.h b/lib/qcd/action/scalar/ScalarInteractionAction.h index 5a322a5e..ca8207bd 100644 --- a/lib/qcd/action/scalar/ScalarInteractionAction.h +++ b/lib/qcd/action/scalar/ScalarInteractionAction.h @@ -37,11 +37,11 @@ directory namespace Grid { // FIXME drop the QCD namespace everywhere here -template -class ScalarInteractionAction : public QCD::Action { -public: + template + class ScalarInteractionAction : public QCD::Action { + public: INHERIT_FIELD_TYPES(Impl); -private: + private: RealD mass_square; RealD lambda; @@ -50,14 +50,19 @@ private: typedef CartesianStencil Stencil; SimpleCompressor compressor; - int npoint = 8; - std::vector directions = {0,1,2,3,0,1,2,3}; // forcing 4 dimensions - std::vector displacements = {1,1,1,1, -1,-1,-1,-1}; + int npoint = 2*Ndim; + std::vector directions;// = {0,1,2,3,0,1,2,3}; // forcing 4 dimensions + std::vector displacements;// = {1,1,1,1, -1,-1,-1,-1}; - public: + public: - ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l){} + ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l), displacements(2*Ndim,0), directions(2*Ndim,0){ + for (int mu = 0 ; mu < Ndim; mu++){ + directions[mu] = mu; directions[mu+Ndim] = mu; + displacements[mu] = 1; displacements[mu+Ndim] = -1; + } + } virtual std::string LogParameters() { std::stringstream sstream; @@ -71,75 +76,74 @@ private: virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} virtual RealD S(const Field &p) { - static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); - phiStencil.HaloExchange(p, compressor); - - Field action(p._grid), pshift(p._grid), phisquared(p._grid); - phisquared = p*p; - action = (2.0*QCD::Nd + mass_square)*phisquared + lambda*phisquared*phisquared; - for (int mu = 0; mu < QCD::Nd; mu++) { - // pshift = Cshift(p, mu, +1); // not efficient, implement with stencils - PARALLEL_FOR_LOOP - for (int i = 0; i < p._grid->oSites(); i++) { - int permute_type; - StencilEntry *SE; - vobj temp2; - vobj *temp; - vobj *t_p; - - SE = phiStencil.GetEntry(permute_type, mu, i); - t_p = &p._odata[i]; - if ( SE->_is_local ) { - temp = &p._odata[SE->_offset]; - if ( SE->_permute ) { - permute(temp2, *temp, permute_type); - action._odata[i] -= temp2*(*t_p) + (*t_p)*temp2; - } else { - action._odata[i] -= *temp*(*t_p) + (*t_p)*(*temp); - } - } else { - action._odata[i] -= phiStencil.CommBuf()[SE->_offset]*(*t_p) + (*t_p)*phiStencil.CommBuf()[SE->_offset]; - } - } - // action -= pshift*p + p*pshift; - } - // NB the trace in the algebra is normalised to 1/2 - // minus sign coming from the antihermitian fields - return -(TensorRemove(sum(trace(action)))).real(); + assert(p._grid->Nd() == Ndim); + static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); + phiStencil.HaloExchange(p, compressor); + Field action(p._grid), pshift(p._grid), phisquared(p._grid); + phisquared = p*p; + action = (2.0*Ndim + mass_square)*phisquared + lambda*phisquared*phisquared; + for (int mu = 0; mu < Ndim; mu++) { + // pshift = Cshift(p, mu, +1); // not efficient, implement with stencils + parallel_for (int i = 0; i < p._grid->oSites(); i++) { + int permute_type; + StencilEntry *SE; + vobj temp2; + vobj *temp; + vobj *t_p; + + SE = phiStencil.GetEntry(permute_type, mu, i); + t_p = &p._odata[i]; + if ( SE->_is_local ) { + temp = &p._odata[SE->_offset]; + if ( SE->_permute ) { + permute(temp2, *temp, permute_type); + action._odata[i] -= temp2*(*t_p) + (*t_p)*temp2; + } else { + action._odata[i] -= *temp*(*t_p) + (*t_p)*(*temp); + } + } else { + action._odata[i] -= phiStencil.CommBuf()[SE->_offset]*(*t_p) + (*t_p)*phiStencil.CommBuf()[SE->_offset]; + } + } + // action -= pshift*p + p*pshift; + } + // NB the trace in the algebra is normalised to 1/2 + // minus sign coming from the antihermitian fields + return -(TensorRemove(sum(trace(action)))).real(); }; virtual void deriv(const Field &p, Field &force) { - force = (2.0*QCD::Nd + mass_square)*p + 2.0*lambda*p*p*p; - // move this outside - static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); - phiStencil.HaloExchange(p, compressor); - - //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); - for (int point = 0; point < npoint; point++) { - PARALLEL_FOR_LOOP - for (int i = 0; i < p._grid->oSites(); i++) { - vobj *temp; - vobj temp2; - int permute_type; - StencilEntry *SE; - SE = phiStencil.GetEntry(permute_type, point, i); - - if ( SE->_is_local ) { - temp = &p._odata[SE->_offset]; - if ( SE->_permute ) { - permute(temp2, *temp, permute_type); - force._odata[i] -= temp2; - } else { - force._odata[i] -= *temp; - } - } else { - force._odata[i] -= phiStencil.CommBuf()[SE->_offset]; - } - } - } + assert(p._grid->Nd() == Ndim); + force = (2.0*Ndim + mass_square)*p + 2.0*lambda*p*p*p; + // move this outside + static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); + phiStencil.HaloExchange(p, compressor); + + //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); + for (int point = 0; point < npoint; point++) { + parallel_for (int i = 0; i < p._grid->oSites(); i++) { + vobj *temp; + vobj temp2; + int permute_type; + StencilEntry *SE; + SE = phiStencil.GetEntry(permute_type, point, i); + + if ( SE->_is_local ) { + temp = &p._odata[SE->_offset]; + if ( SE->_permute ) { + permute(temp2, *temp, permute_type); + force._odata[i] -= temp2; + } else { + force._odata[i] -= *temp; + } + } else { + force._odata[i] -= phiStencil.CommBuf()[SE->_offset]; + } + } + } } -}; - + }; + } // namespace Grid #endif // SCALAR_INT_ACTION_H diff --git a/lib/qcd/hmc/GenericHMCrunner.h b/lib/qcd/hmc/GenericHMCrunner.h index 353b4905..4f6c1af0 100644 --- a/lib/qcd/hmc/GenericHMCrunner.h +++ b/lib/qcd/hmc/GenericHMCrunner.h @@ -210,6 +210,9 @@ typedef HMCWrapperTemplate typedef HMCWrapperTemplate ScalarAdjGenericHMCRunner; +template +using ScalarNxNAdjGenericHMCRunner = HMCWrapperTemplate < ScalarNxNAdjImplR, MinimumNorm2, ScalarNxNMatrixFields >; + } // namespace QCD } // namespace Grid diff --git a/lib/qcd/representations/hmc_types.h b/lib/qcd/representations/hmc_types.h index b4991941..3fee377e 100644 --- a/lib/qcd/representations/hmc_types.h +++ b/lib/qcd/representations/hmc_types.h @@ -64,6 +64,9 @@ typedef Representations NoHirep; typedef Representations > ScalarFields; typedef Representations > ScalarMatrixFields; +template < int Colours> +using ScalarNxNMatrixFields = Representations::Field> >; + // Helper classes to access the elements // Strips the first N parameters from the tuple // sequence of classes to obtain the S sequence diff --git a/lib/stencil/Stencil.h b/lib/stencil/Stencil.h index d1c28e78..887142c4 100644 --- a/lib/stencil/Stencil.h +++ b/lib/stencil/Stencil.h @@ -286,7 +286,7 @@ class CartesianStencil { // Stencil runs along coordinate axes only; NO diagonal { int dimension = _directions[point]; int displacement = _distances[point]; - + int fd = _grid->_fdimensions[dimension]; int rd = _grid->_rdimensions[dimension]; diff --git a/tests/hmc/Test_hmc_ScalarActionNxN.cc b/tests/hmc/Test_hmc_ScalarActionNxN.cc index f63936b5..b3ce6840 100644 --- a/tests/hmc/Test_hmc_ScalarActionNxN.cc +++ b/tests/hmc/Test_hmc_ScalarActionNxN.cc @@ -32,68 +32,116 @@ class ScalarActionParameters : Serializable { GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarActionParameters, double, mass_squared, double, lambda); + + template + ScalarActionParameters(Reader& Reader){ + read(Reader, "ScalarAction", *this); + } + }; } int main(int argc, char **argv) { using namespace Grid; using namespace Grid::QCD; - + typedef Grid::JSONReader Serialiser; + Grid_init(&argc, &argv); int threads = GridThread::GetThreads(); // here make a routine to print all the relevant information on the run std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; // Typedefs to simplify notation - typedef ScalarAdjGenericHMCRunner HMCWrapper; // Uses the default minimum norm, real scalar fields - + constexpr int Ncolours = 4; + constexpr int Ndimensions = 3; + typedef ScalarNxNAdjGenericHMCRunner HMCWrapper; // Uses the default minimum norm, real scalar fields + typedef ScalarAdjActionR ScalarAction; //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: HMCWrapper TheHMC; + TheHMC.ReadCommandLine(argc, argv); + + if (TheHMC.ParameterFile.empty()){ + std::cout << "Input file not specified." + << "Use --ParameterFile option in the command line.\nAborting" + << std::endl; + exit(1); + } + Serialiser Reader(TheHMC.ParameterFile); // Grid from the command line GridModule ScalarGrid; - ScalarGrid.set_full(SpaceTimeGrid::makeFourDimGrid( - GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), - GridDefaultMpi())); - ScalarGrid.set_rb(SpaceTimeGrid::makeFourDimRedBlackGrid(ScalarGrid.get_full())); + if (GridDefaultLatt().size() != Ndimensions){ + std::cout << "Incorrect dimension of the grid\n. Expected dim="<< Ndimensions << std::endl; + exit(1); + } + if (GridDefaultMpi().size() != Ndimensions){ + std::cout << "Incorrect dimension of the mpi grid\n. Expected dim="<< Ndimensions << std::endl; + exit(1); + } + ScalarGrid.set_full(new GridCartesian(GridDefaultLatt(),GridDefaultSimd(Ndimensions, vComplex::Nsimd()),GridDefaultMpi())); + ScalarGrid.set_rb(new GridRedBlackCartesian(ScalarGrid.get_full())); TheHMC.Resources.AddGrid("scalar", ScalarGrid); - // Possibile to create the module by hand - // hardcoding parameters or using a Reader + std::cout << "Lattice size : " << GridDefaultLatt() << std::endl; // Checkpointer definition - CheckpointerParameters CPparams; - CPparams.config_prefix = "ckpoint_scalar_lat"; - CPparams.rng_prefix = "ckpoint_scalar_rng"; - CPparams.saveInterval = 50; - CPparams.format = "IEEE64BIG"; - + CheckpointerParameters CPparams(Reader); TheHMC.Resources.LoadBinaryCheckpointer(CPparams); - RNGModuleParameters RNGpar; - RNGpar.serial_seeds = "1 2 3 4 5"; - RNGpar.parallel_seeds = "6 7 8 9 10"; + RNGModuleParameters RNGpar(Reader); TheHMC.Resources.SetRNGSeeds(RNGpar); ///////////////////////////////////////////////////////////// // Collect actions, here use more encapsulation // Scalar action in adjoint representation - ScalarActionParameters SPar; - SPar.mass_squared = 0.5; - SPar.lambda = 0.1; - ScalarAdjActionR Saction(SPar.mass_squared, SPar.lambda); + ScalarActionParameters SPar(Reader); + ScalarAction Saction(SPar.mass_squared, SPar.lambda); // Collect actions - ActionLevel Level1(1); + ActionLevel> Level1(1); Level1.push_back(&Saction); TheHMC.TheAction.push_back(Level1); ///////////////////////////////////////////////////////////// + TheHMC.Parameters.initialize(Reader); - // HMC parameters are serialisable - TheHMC.Parameters.MD.MDsteps = 20; - TheHMC.Parameters.MD.trajL = 1.0; - - TheHMC.ReadCommandLine(argc, argv); TheHMC.Run(); Grid_finalize(); } // main + +/* Examples for input files + +JSON + +{ + "Checkpointer": { + "config_prefix": "ckpoint_scalar_lat", + "rng_prefix": "ckpoint_scalar_rng", + "saveInterval": 1, + "format": "IEEE64BIG" + }, + "RandomNumberGenerator": { + "serial_seeds": "1 2 3 4 6", + "parallel_seeds": "6 7 8 9 11" + }, + "ScalarAction":{ + "mass_squared": 0.5, + "lambda": 0.1 + }, + "HMC":{ + "StartTrajectory": 0, + "Trajectories": 100, + "MetropolisTest": true, + "NoMetropolisUntil": 10, + "StartingType": "HotStart", + "MD":{ + "name": "MinimumNorm2", + "MDsteps": 15, + "trajL": 2.0 + } + } +} + + +XML example not provided yet + +*/ \ No newline at end of file From 43c817cc67c6447bbf69bfc7d7772fba4e7ff9eb Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 11 May 2017 00:07:17 +0100 Subject: [PATCH 6/8] Scalar action: const fix --- lib/qcd/action/scalar/ScalarInteractionAction.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/lib/qcd/action/scalar/ScalarInteractionAction.h b/lib/qcd/action/scalar/ScalarInteractionAction.h index ca8207bd..5f4c630c 100644 --- a/lib/qcd/action/scalar/ScalarInteractionAction.h +++ b/lib/qcd/action/scalar/ScalarInteractionAction.h @@ -88,8 +88,7 @@ namespace Grid { int permute_type; StencilEntry *SE; vobj temp2; - vobj *temp; - vobj *t_p; + const vobj *temp, *t_p; SE = phiStencil.GetEntry(permute_type, mu, i); t_p = &p._odata[i]; @@ -122,7 +121,7 @@ namespace Grid { //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); for (int point = 0; point < npoint; point++) { parallel_for (int i = 0; i < p._grid->oSites(); i++) { - vobj *temp; + const vobj *temp; vobj temp2; int permute_type; StencilEntry *SE; From d1ece741370d1b829f5946afc7c21c585a158d31 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 11 May 2017 11:40:44 +0100 Subject: [PATCH 7/8] HMC scalar test: magnetisation measurement --- tests/hmc/Test_hmc_ScalarActionNxN.cc | 54 ++++++++++++++++++++++++--- 1 file changed, 49 insertions(+), 5 deletions(-) diff --git a/tests/hmc/Test_hmc_ScalarActionNxN.cc b/tests/hmc/Test_hmc_ScalarActionNxN.cc index b3ce6840..bcaee31d 100644 --- a/tests/hmc/Test_hmc_ScalarActionNxN.cc +++ b/tests/hmc/Test_hmc_ScalarActionNxN.cc @@ -39,11 +39,50 @@ class ScalarActionParameters : Serializable { } }; - } + +using namespace Grid; +using namespace Grid::QCD; + +template +class MagLogger : public HmcObservable { +public: + typedef typename Impl::Field Field; + typedef typename Impl::Simd::scalar_type Trace; + + void TrajectoryComplete(int traj, + Field &U, + GridSerialRNG &sRNG, + GridParallelRNG &pRNG) { + + int def_prec = std::cout.precision(); + + std::cout << std::setprecision(std::numeric_limits::digits10 + 1); + std::cout << GridLogMessage + << "m= " << TensorRemove(trace(sum(U))) << std::endl; + std::cout << GridLogMessage + << "m^2= " << TensorRemove(trace(sum(U)*sum(U))) << std::endl; + std::cout.precision(def_prec); + + } +private: + +}; + +template +class MagMod: public ObservableModule, NoParameters>{ + typedef ObservableModule, NoParameters> ObsBase; + using ObsBase::ObsBase; // for constructors + + // acquire resource + virtual void initialize(){ + this->ObservablePtr.reset(new MagLogger()); + } +public: + MagMod(): ObsBase(NoParameters()){} +}; + int main(int argc, char **argv) { - using namespace Grid; - using namespace Grid::QCD; typedef Grid::JSONReader Serialiser; Grid_init(&argc, &argv); @@ -52,7 +91,7 @@ int main(int argc, char **argv) { std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; // Typedefs to simplify notation - constexpr int Ncolours = 4; + constexpr int Ncolours = 2; constexpr int Ndimensions = 3; typedef ScalarNxNAdjGenericHMCRunner HMCWrapper; // Uses the default minimum norm, real scalar fields typedef ScalarAdjActionR ScalarAction; @@ -89,6 +128,11 @@ int main(int argc, char **argv) { RNGModuleParameters RNGpar(Reader); TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + typedef MagMod MagObs; + TheHMC.Resources.AddObservable(); + ///////////////////////////////////////////////////////////// // Collect actions, here use more encapsulation @@ -144,4 +188,4 @@ JSON XML example not provided yet -*/ \ No newline at end of file +*/ From 3f858d675557536feb6bac6312e4205c987857d9 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 May 2017 13:25:14 +0200 Subject: [PATCH 8/8] Scalar: phi^2 observable --- tests/hmc/Test_hmc_ScalarActionNxN.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/hmc/Test_hmc_ScalarActionNxN.cc b/tests/hmc/Test_hmc_ScalarActionNxN.cc index bcaee31d..a7490f51 100644 --- a/tests/hmc/Test_hmc_ScalarActionNxN.cc +++ b/tests/hmc/Test_hmc_ScalarActionNxN.cc @@ -62,6 +62,8 @@ public: << "m= " << TensorRemove(trace(sum(U))) << std::endl; std::cout << GridLogMessage << "m^2= " << TensorRemove(trace(sum(U)*sum(U))) << std::endl; + std::cout << GridLogMessage + << "phi^2= " << TensorRemove(sum(trace(U*U))) << std::endl; std::cout.precision(def_prec); }