From b32fd473f819516f9193eac130548f86e58bd1d2 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:28:28 +0200 Subject: [PATCH] Start at the Domain decomposed supprt --- .../pseudofermion/DomainDecomposedBoundary.h | 346 ++++++++++++++++++ 1 file changed, 346 insertions(+) create mode 100644 Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h diff --git a/Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h b/Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h new file mode 100644 index 00000000..36a849e0 --- /dev/null +++ b/Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h @@ -0,0 +1,346 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: Peter Boyle +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 + 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 */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +/////////////////////////////////////// +// Two flavour ratio +/////////////////////////////////////// +template +class DomainDecomposedBoundary { +public: + INHERIT_IMPL_TYPES(Impl); + + typedef typename GaugeField::vector_type vector_type; //SIMD-vectorized complex type + typedef typename GaugeField::scalar_type scalar_type; //scalar complex type + + typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction + typedef iScalar > > ScalarType; //complex phase for each site + typedef Lattice LatticeLorentzScalarType; + typedef Lattice LatticeScalarType; + + Coordinate Block; + DDHMCFilter Filter; + const int Omega=0; + const int OmegaBar=1; + + void ProjectBoundaryBothDomains (FermionField &f,int sgn) + { + assert((sgn==1)||(sgn==-1)); + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + + GridBase *grid = f.Grid(); + LatticeInteger coor(grid); + LatticeInteger face(grid); + LatticeInteger nface(grid); nface=Zero(); + + ComplexField zz(grid); zz=Zero(); + + FermionField projected(grid); projected=Zero(); + FermionField sp_proj (grid); + + int dims = grid->Nd(); + int isDWF= (dims==Nd+1); + assert((dims==Nd)||(dims==Nd+1)); + + for(int mu=0;mu1,f,projected); + } + void ProjectDomain(FermionField &f,int cb) + { + GridBase *grid = f.Grid(); + ComplexField zz(grid); zz=Zero(); + LatticeInteger coor(grid); + LatticeInteger domaincb(grid); domaincb=Zero(); + for(int d=0;dNd();d++){ + LatticeCoordinate(coor,mu); + domaincb = domaincb + div(coor,Block[d]); + } + f = where(mod(domaincb,2)==cb,f,zz); + }; + + void ProjectOmegaBar (FermionField &f) {ProjectDomain(f,OmegaBar);} + void ProjectOmega (FermionField &f) {ProjectDomain(f,Omega);} + + // See my notes(!). + // Notation: Following Luscher, we introduce projectors $\hPdb$ with both spinor and space structure + // projecting all spinor elements in $\Omega$ connected by $\Ddb$ to $\bar{\Omega}$, + void ProjectBoundaryBar(FermionField &f) + { + ProjectBoundaryBothDomains(f); + ProjectOmega(f); + } + // and $\hPd$ projecting all spinor elements in $\bar{\Omega}$ connected by $\Dd$ to $\Omega$. + void ProjectBoundary (FermionField &f) + { + ProjectBoundaryBothDomains(f); + ProjectOmegaBar(f); + }; + + void dBoundary (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmegaBar(tmp); + Op.M(tmp,out); + ProjectOmega(out); + }; + void dBoundaryBar (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmega(tmp); + Op.M(tmp,out); + ProjectOmegaBar(out); + }; + void dOmega (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmega(tmp); + Op.M(tmp,out); + ProjectOmega(out); + }; + void dOmegaBar (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmegaBar(tmp); + Op.M(tmp,out); + ProjectOmegaBar(out); + }; + + void SolveOmega (FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + void SolveOmegaBar(FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + void SolveOmegaAndOmegaBar(FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + void dInverse (FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + + // R = Pdbar - Pdbar DomegaInv Dd DomegabarInv Ddbar + void R(FermionOperator &Op,FermionOperator &OpDirichlet,FermionField &in,FermionField &out) + { + FermionField tmp1(Op.FermionGrid()); + FermionField tmp2(Op.FermionGrid()); + dBoundaryBar(Op,in,tmp1); + SolveOmegaBar(OpDirichlet,tmp1,tmp2); // 1/2 cost + dBoundary(Op,tmp2,tmp1); + SolveOmega(OpDirichlet,tmp1,tmp2); // 1/2 cost + out = in - tmp2 ; + ProjectBoundaryBar(out); + }; + + // R = Pdbar - Pdbar Dinv Ddbar + void Rinverse(FermionField &in,FermionField &out) + { + FermionField tmp1(NumOp.FermionGrid()); + out = in; + ProjectBoundaryBar(out); + dInverse(out,tmp1); + ProjectBoundaryBar(tmp1); + out = out -tmp1; + }; + +} + +template +class DomainDecomposedBoundaryPseudoFermionAction : public Action { +public: + INHERIT_IMPL_TYPES(Impl); + +private: + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + FermionOperator & NumOpDirichlet;// the basic operator + FermionOperator & DenOpDirichlet;// the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; + + FermionField Phi; // the pseudo fermion field for this trajectory + + +public: + DomainBoundaryPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + FermionOperator &_NumOpDirichlet, + FermionOperator &_DenOpDirichlet, + OperatorFunction & DS, + OperatorFunction & AS, + Coordinate &_Block + ) : NumOp(_NumOp), + DenOp(_DenOp), + DerivativeSolver(DS), + ActionSolver(AS), + Phi(_NumOp.FermionGrid()), + Block(_Block) + // LinkFilter(Block) + {}; + + virtual std::string action_name(){return "DomainBoundaryPseudoFermionRatioAction";} + + virtual std::string LogParameters(){ + std::stringstream sstream; + sstream << GridLogMessage << "["< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707.... + // + RealD scale = std::sqrt(0.5); + + FermionField eta(NumOp.FermionGrid()); + FermionField tmp(NumOp.FermionGrid()); + + gaussian(pRNG,eta); + + ProjectBoundary(eta); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + // Note: this hard codes normal equations type solvers; alternate implementation needed for + // non-herm style solvers. + MdagMLinearOperator ,FermionField> MdagMOp(NumOp); + + DenOp.Mdag(eta,Phi); // Mdag eta + tmp = Zero(); + ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta + NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta + + Phi=Phi*scale; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag V (Mdag M)^-1 Vdag phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + NumOp.Mdag(Phi,Y); // Y= Vdag phi + X=Zero(); + ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi + + RealD action = norm2(Y); + + return action; + }; + + ////////////////////////////////////////////////////// + // dS/du = phi^dag dV (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi + // + phi^dag V (Mdag M)^-1 dV^dag phi + ////////////////////////////////////////////////////// + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + GaugeField force(NumOp.GaugeGrid()); + + + //Y=Vdag phi + //X = (Mdag M)^-1 V^dag phi + //Y = (Mdag)^-1 V^dag phi + NumOp.Mdag(Phi,Y); // Y= Vdag phi + X=Zero(); + DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi + + // phi^dag V (Mdag M)^-1 dV^dag phi + NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; + + // phi^dag dV (Mdag M)^-1 V^dag phi + NumOp.MDeriv(force , Phi, X ,DaggerNo ); dSdU=dSdU+force; + + // - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi + DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force; + DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; + + dSdU *= -1.0; + //dSdU = - Ta(dSdU); + + }; +}; + +NAMESPACE_END(Grid); + +