portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-11 14:40:46 +01:00
Code Releases Activity

Namespace QCD gone

Browse Source
This commit is contained in:
paboyle 2018-01-14 22:38:22 +00:00
parent cc349c6512
commit 6ee5ea6b32
1 changed files with 84 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

170
lib/qcd/action/scalar/ScalarInteractionAction.h
View File

@ -24,125 +24,123 @@
See the full license in the file "LICENSE" in the top level distribution See the full license in the file "LICENSE" in the top level distribution
directory directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#ifndef SCALAR_INT_ACTION_H #ifndef SCALAR_INT_ACTION_H
#define SCALAR_INT_ACTION_H #define SCALAR_INT_ACTION_H
// Note: this action can completely absorb the ScalarAction for real float fields // Note: this action can completely absorb the ScalarAction for real float fields
// use the scalarObjs to generalise the structure // use the scalarObjs to generalise the structure
namespace Grid { NAMESPACE_BEGIN(Grid);
// FIXME drop the QCD namespace everywhere here
template <class Impl, int Ndim > template <class Impl, int Ndim >
class ScalarInteractionAction : public QCD::Action<typename Impl::Field> { class ScalarInteractionAction : public QCD::Action<typename Impl::Field> {
public: public:
INHERIT_FIELD_TYPES(Impl); INHERIT_FIELD_TYPES(Impl);
private: private:
RealD mass_square; RealD mass_square;
RealD lambda; RealD lambda;
typedef typename Field::vector_object vobj; typedef typename Field::vector_object vobj;
typedef CartesianStencil<vobj,vobj> Stencil; typedef CartesianStencil<vobj,vobj> Stencil;
SimpleCompressor<vobj> compressor; SimpleCompressor<vobj> compressor;
int npoint = 2*Ndim; int npoint = 2*Ndim;
std::vector<int> directions;// = {0,1,2,3,0,1,2,3}; // forcing 4 dimensions std::vector<int> directions;// = {0,1,2,3,0,1,2,3}; // forcing 4 dimensions
std::vector<int> displacements;// = {1,1,1,1, -1,-1,-1,-1}; std::vector<int> displacements;// = {1,1,1,1, -1,-1,-1,-1};
public: public:
ScalarInteractionAction(RealD ms, RealD l) : mass_square(ms), lambda(l), displacements(2*Ndim,0), directions(2*Ndim,0){ 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++){ for (int mu = 0 ; mu < Ndim; mu++){
directions[mu] = mu; directions[mu+Ndim] = mu; directions[mu] = mu; directions[mu+Ndim] = mu;
displacements[mu] = 1; displacements[mu+Ndim] = -1; displacements[mu] = 1; displacements[mu+Ndim] = -1;
}
} }
}
virtual std::string LogParameters() { virtual std::string LogParameters() {
std::stringstream sstream; std::stringstream sstream;
sstream << GridLogMessage << "[ScalarAction] lambda : " << lambda << std::endl; sstream << GridLogMessage << "[ScalarAction] lambda : " << lambda << std::endl;
sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl; sstream << GridLogMessage << "[ScalarAction] mass_square : " << mass_square << std::endl;
return sstream.str(); return sstream.str();
} }
virtual std::string action_name() {return "ScalarAction";} virtual std::string action_name() {return "ScalarAction";}
virtual void refresh(const Field &U, GridParallelRNG &pRNG) {} virtual void refresh(const Field &U, GridParallelRNG &pRNG) {}
virtual RealD S(const Field &p) { virtual RealD S(const Field &p) {
assert(p._grid->Nd() == Ndim); assert(p._grid->Nd() == Ndim);
static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); static Stencil phiStencil(p._grid, npoint, 0, directions, displacements);
phiStencil.HaloExchange(p, compressor); phiStencil.HaloExchange(p, compressor);
Field action(p._grid), pshift(p._grid), phisquared(p._grid); Field action(p._grid), pshift(p._grid), phisquared(p._grid);
phisquared = p*p; phisquared = p*p;
action = (2.0*Ndim + mass_square)*phisquared - lambda/24.*phisquared*phisquared; action = (2.0*Ndim + mass_square)*phisquared - lambda/24.*phisquared*phisquared;
for (int mu = 0; mu < Ndim; mu++) { for (int mu = 0; mu < Ndim; mu++) {
// pshift = Cshift(p, mu, +1); // not efficient, implement with stencils // pshift = Cshift(p, mu, +1); // not efficient, implement with stencils
parallel_for (int i = 0; i < p._grid->oSites(); i++) { parallel_for (int i = 0; i < p._grid->oSites(); i++) {
int permute_type; int permute_type;
StencilEntry *SE; StencilEntry *SE;
vobj temp2; vobj temp2;
const vobj *temp, *t_p; const vobj *temp, *t_p;
SE = phiStencil.GetEntry(permute_type, mu, i); SE = phiStencil.GetEntry(permute_type, mu, i);
t_p = &p._odata[i]; t_p = &p._odata[i];
if ( SE->_is_local ) { if ( SE->_is_local ) {
temp = &p._odata[SE->_offset]; temp = &p._odata[SE->_offset];
if ( SE->_permute ) { if ( SE->_permute ) {
permute(temp2, *temp, permute_type); permute(temp2, *temp, permute_type);
action._odata[i] -= temp2*(*t_p) + (*t_p)*temp2; action._odata[i] -= temp2*(*t_p) + (*t_p)*temp2;
} else {
action._odata[i] -= (*temp)*(*t_p) + (*t_p)*(*temp);
}
} else { } else {
action._odata[i] -= phiStencil.CommBuf()[SE->_offset]*(*t_p) + (*t_p)*phiStencil.CommBuf()[SE->_offset]; 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 // action -= pshift*p + p*pshift;
// minus sign coming from the antihermitian fields }
return -(TensorRemove(sum(trace(action)))).real(); // 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) {
assert(p._grid->Nd() == Ndim); assert(p._grid->Nd() == Ndim);
force = (2.0*Ndim + mass_square)*p - lambda/12.*p*p*p; force = (2.0*Ndim + mass_square)*p - lambda/12.*p*p*p;
// move this outside // move this outside
static Stencil phiStencil(p._grid, npoint, 0, directions, displacements); static Stencil phiStencil(p._grid, npoint, 0, directions, displacements);
phiStencil.HaloExchange(p, compressor); phiStencil.HaloExchange(p, compressor);
//for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1); //for (int mu = 0; mu < QCD::Nd; mu++) force -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
for (int point = 0; point < npoint; point++) { for (int point = 0; point < npoint; point++) {
parallel_for (int i = 0; i < p._grid->oSites(); i++) { parallel_for (int i = 0; i < p._grid->oSites(); i++) {
const vobj *temp; const vobj *temp;
vobj temp2; vobj temp2;
int permute_type; int permute_type;
StencilEntry *SE; StencilEntry *SE;
SE = phiStencil.GetEntry(permute_type, point, i); SE = phiStencil.GetEntry(permute_type, point, i);
if ( SE->_is_local ) { if ( SE->_is_local ) {
temp = &p._odata[SE->_offset]; temp = &p._odata[SE->_offset];
if ( SE->_permute ) { if ( SE->_permute ) {
permute(temp2, *temp, permute_type); permute(temp2, *temp, permute_type);
force._odata[i] -= temp2; force._odata[i] -= temp2;
} else {
force._odata[i] -= *temp;
}
} else { } else {
force._odata[i] -= phiStencil.CommBuf()[SE->_offset]; force._odata[i] -= *temp;
} }
} else {
force._odata[i] -= phiStencil.CommBuf()[SE->_offset];
} }
} }
} }
}; }
};
} // namespace Grid NAMESPACE_END(Grid);
#endif // SCALAR_INT_ACTION_H #endif // SCALAR_INT_ACTION_H