mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-07 09:57:07 +01:00
Compare commits
41 Commits
e652fc2825
...
21de6f7da8
Author | SHA1 | Date | |
---|---|---|---|
|
21de6f7da8 | ||
|
dbe39f9ce0 | ||
|
ab3de50d5e | ||
|
c545bd2139 | ||
|
6a1c64fbdd | ||
|
b75809ed61 | ||
|
ecaf228e5c | ||
|
6d015ae8fc | ||
|
233150d93f | ||
|
7af8c77a52 | ||
|
a957e7bfa1 | ||
|
cee4c8ce8c | ||
|
96bf814d8c | ||
|
7ddc422788 | ||
|
e465fce201 | ||
|
d41542c64b | ||
|
785bc7a14f | ||
|
1a1fe85428 | ||
|
0000d2e558 | ||
|
b1ba209696 | ||
|
cb3e529b1e | ||
|
717f647418 | ||
|
98e7418187 | ||
|
fe05bf48b1 | ||
|
d2dd8f54e2 | ||
|
7726ee4b16 | ||
|
8729c46169 | ||
|
09f81fe7c3 | ||
|
1876e5b7c0 | ||
|
355ec76257 | ||
|
4f17c8d081 | ||
|
aaab753982 | ||
|
570b72a47b | ||
|
a5798a89ed | ||
|
f7e2f9a401 | ||
|
2848a9b558 | ||
|
d4868991af | ||
|
e99d42404e | ||
|
3ba019c747 | ||
|
47429218bb | ||
|
8d305df0db |
@ -245,9 +245,10 @@ until convergence
|
||||
_HermOp(src_n,tmp);
|
||||
// std::cout << GridLogMessage<< tmp<<std::endl; exit(0);
|
||||
// std::cout << GridLogIRL << " _HermOp " << norm2(tmp) << std::endl;
|
||||
RealD vnum = real(innerProduct(src_n,tmp)); // HermOp.
|
||||
// RealD vnum = real(innerProduct(src_n,tmp)); // HermOp.
|
||||
RealD vnum = real(innerProduct(tmp,tmp)); // HermOp^2.
|
||||
RealD vden = norm2(src_n);
|
||||
RealD na = vnum/vden;
|
||||
RealD na = std::sqrt(vnum/vden);
|
||||
if (fabs(evalMaxApprox/na - 1.0) < 0.0001)
|
||||
i=_MAX_ITER_IRL_MEVAPP_;
|
||||
evalMaxApprox = na;
|
||||
@ -255,6 +256,7 @@ until convergence
|
||||
src_n = tmp;
|
||||
}
|
||||
}
|
||||
std::cout << GridLogIRL << " Final evalMaxApprox " << evalMaxApprox << std::endl;
|
||||
|
||||
std::vector<RealD> lme(Nm);
|
||||
std::vector<RealD> lme2(Nm);
|
||||
|
@ -510,7 +510,6 @@ public:
|
||||
grid->SendToRecvFromBegin(fwd_req,
|
||||
(void *)&hsend_buf[d*buffer_size], xmit_to_rank,
|
||||
(void *)&hrecv_buf[d*buffer_size], recv_from_rank, bytes, tag);
|
||||
acceleratorCopyToDevice(&hrecv_buf[d*buffer_size],&recv_buf[d*buffer_size],bytes);
|
||||
#endif
|
||||
t_comms+=usecond()-t;
|
||||
}
|
||||
@ -531,7 +530,6 @@ public:
|
||||
grid->SendToRecvFromBegin(bwd_req,
|
||||
(void *)&hsend_buf[(d+depth)*buffer_size], recv_from_rank,
|
||||
(void *)&hrecv_buf[(d+depth)*buffer_size], xmit_to_rank, bytes,tag);
|
||||
acceleratorCopyToDevice(&hrecv_buf[(d+depth)*buffer_size],&recv_buf[(d+depth)*buffer_size],bytes);
|
||||
#endif
|
||||
t_comms+=usecond()-t;
|
||||
}
|
||||
@ -555,8 +553,13 @@ public:
|
||||
|
||||
t=usecond();
|
||||
grid->CommsComplete(fwd_req);
|
||||
#ifndef ACCELERATOR_AWARE_MPI
|
||||
for ( int d=0;d < depth ; d ++ ) {
|
||||
acceleratorCopyToDevice(&hrecv_buf[d*buffer_size],&recv_buf[d*buffer_size],bytes);
|
||||
}
|
||||
#endif
|
||||
t_comms+= usecond() - t;
|
||||
|
||||
|
||||
t=usecond();
|
||||
for ( int d=0;d < depth ; d ++ ) {
|
||||
ScatterSlice(recv_buf,to,nld-depth+d,dimension,plane*buffer_size); plane++;
|
||||
@ -565,6 +568,11 @@ public:
|
||||
|
||||
t=usecond();
|
||||
grid->CommsComplete(bwd_req);
|
||||
#ifndef ACCELERATOR_AWARE_MPI
|
||||
for ( int d=0;d < depth ; d ++ ) {
|
||||
acceleratorCopyToDevice(&hrecv_buf[(d+depth)*buffer_size],&recv_buf[(d+depth)*buffer_size],bytes);
|
||||
}
|
||||
#endif
|
||||
t_comms+= usecond() - t;
|
||||
|
||||
t=usecond();
|
||||
|
196
Grid/qcd/action/fermion/CompactWilsonCloverFermion5D.h
Normal file
196
Grid/qcd/action/fermion/CompactWilsonCloverFermion5D.h
Normal file
@ -0,0 +1,196 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/CompactWilsonCloverFermion5D.h
|
||||
|
||||
Copyright (C) 2020 - 2025
|
||||
|
||||
Author: Daniel Richtmann <daniel.richtmann@gmail.com>
|
||||
Author: Nils Meyer <nils.meyer@ur.de>
|
||||
Author: Christoph Lehner <christoph@lhnr.de>
|
||||
|
||||
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
|
||||
|
||||
#include <Grid/qcd/action/fermion/WilsonFermion5D.h>
|
||||
#include <Grid/qcd/action/fermion/WilsonCloverTypes.h>
|
||||
#include <Grid/qcd/action/fermion/WilsonCloverHelpers.h>
|
||||
#include <Grid/qcd/action/fermion/CloverHelpers.h>
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
// see Grid/qcd/action/fermion/CompactWilsonCloverFermion.h for description
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
class CompactWilsonCloverFermion5D : public WilsonFermion5D<Impl>,
|
||||
public WilsonCloverHelpers<Impl>,
|
||||
public CompactWilsonCloverHelpers<Impl> {
|
||||
/////////////////////////////////////////////
|
||||
// Sizes
|
||||
/////////////////////////////////////////////
|
||||
|
||||
public:
|
||||
|
||||
INHERIT_COMPACT_CLOVER_SIZES(Impl);
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Type definitions
|
||||
/////////////////////////////////////////////
|
||||
|
||||
public:
|
||||
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
INHERIT_CLOVER_TYPES(Impl);
|
||||
INHERIT_COMPACT_CLOVER_TYPES(Impl);
|
||||
|
||||
typedef WilsonFermion5D<Impl> WilsonBase;
|
||||
typedef WilsonCloverHelpers<Impl> Helpers;
|
||||
typedef CompactWilsonCloverHelpers<Impl> CompactHelpers;
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Constructors
|
||||
/////////////////////////////////////////////
|
||||
|
||||
public:
|
||||
|
||||
CompactWilsonCloverFermion5D(GaugeField& _Umu,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
const RealD _mass,
|
||||
const RealD _csw_r = 0.0,
|
||||
const RealD _csw_t = 0.0,
|
||||
const RealD _cF = 1.0,
|
||||
const ImplParams& impl_p = ImplParams());
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Member functions (implementing interface)
|
||||
/////////////////////////////////////////////
|
||||
|
||||
public:
|
||||
|
||||
virtual void Instantiatable() {};
|
||||
int ConstEE() override { return 0; };
|
||||
int isTrivialEE() override { return 0; };
|
||||
|
||||
void Dhop(const FermionField& in, FermionField& out, int dag) override;
|
||||
|
||||
void DhopOE(const FermionField& in, FermionField& out, int dag) override;
|
||||
|
||||
void DhopEO(const FermionField& in, FermionField& out, int dag) override;
|
||||
|
||||
void DhopDir(const FermionField& in, FermionField& out, int dir, int disp) override;
|
||||
|
||||
void DhopDirAll(const FermionField& in, std::vector<FermionField>& out) /* override */;
|
||||
|
||||
void M(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void Mdag(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void Meooe(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void MeooeDag(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void Mooee(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void MooeeDag(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void MooeeInv(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void MooeeInvDag(const FermionField& in, FermionField& out) override;
|
||||
|
||||
void Mdir(const FermionField& in, FermionField& out, int dir, int disp) override;
|
||||
|
||||
void MdirAll(const FermionField& in, std::vector<FermionField>& out) override;
|
||||
|
||||
void MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) override;
|
||||
|
||||
void MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) override;
|
||||
|
||||
void MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) override;
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Member functions (internals)
|
||||
/////////////////////////////////////////////
|
||||
|
||||
void MooeeInternal(const FermionField& in,
|
||||
FermionField& out,
|
||||
const CloverDiagonalField& diagonal,
|
||||
const CloverTriangleField& triangle);
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Helpers
|
||||
/////////////////////////////////////////////
|
||||
|
||||
void ImportGauge(const GaugeField& _Umu) override;
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Helpers
|
||||
/////////////////////////////////////////////
|
||||
|
||||
private:
|
||||
|
||||
template<class Field>
|
||||
const MaskField* getCorrectMaskField(const Field &in) const {
|
||||
if(in.Grid()->_isCheckerBoarded) {
|
||||
if(in.Checkerboard() == Odd) {
|
||||
return &this->BoundaryMaskOdd;
|
||||
} else {
|
||||
return &this->BoundaryMaskEven;
|
||||
}
|
||||
} else {
|
||||
return &this->BoundaryMask;
|
||||
}
|
||||
}
|
||||
|
||||
template<class Field>
|
||||
void ApplyBoundaryMask(Field& f) {
|
||||
const MaskField* m = getCorrectMaskField(f); assert(m != nullptr);
|
||||
assert(m != nullptr);
|
||||
CompactHelpers::ApplyBoundaryMask(f, *m);
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////
|
||||
// Member Data
|
||||
/////////////////////////////////////////////
|
||||
|
||||
public:
|
||||
|
||||
RealD csw_r;
|
||||
RealD csw_t;
|
||||
RealD cF;
|
||||
int n_rhs;
|
||||
|
||||
bool fixedBoundaries;
|
||||
|
||||
CloverDiagonalField Diagonal, DiagonalEven, DiagonalOdd;
|
||||
CloverDiagonalField DiagonalInv, DiagonalInvEven, DiagonalInvOdd;
|
||||
|
||||
CloverTriangleField Triangle, TriangleEven, TriangleOdd;
|
||||
CloverTriangleField TriangleInv, TriangleInvEven, TriangleInvOdd;
|
||||
|
||||
FermionField Tmp;
|
||||
|
||||
MaskField BoundaryMask, BoundaryMaskEven, BoundaryMaskOdd;
|
||||
};
|
||||
|
||||
NAMESPACE_END(Grid);
|
@ -55,6 +55,7 @@ NAMESPACE_CHECK(Wilson);
|
||||
NAMESPACE_CHECK(WilsonTM);
|
||||
#include <Grid/qcd/action/fermion/WilsonCloverFermion.h> // 4d wilson clover fermions
|
||||
#include <Grid/qcd/action/fermion/CompactWilsonCloverFermion.h> // 4d compact wilson clover fermions
|
||||
#include <Grid/qcd/action/fermion/CompactWilsonCloverFermion5D.h> // 5d compact wilson clover fermions
|
||||
NAMESPACE_CHECK(WilsonClover);
|
||||
#include <Grid/qcd/action/fermion/WilsonFermion5D.h> // 5d base used by all 5d overlap types
|
||||
NAMESPACE_CHECK(Wilson5D);
|
||||
@ -164,12 +165,17 @@ typedef WilsonClover<WilsonTwoIndexAntiSymmetricImplD> WilsonCloverTwoIndexAntiS
|
||||
|
||||
// Compact Clover fermions
|
||||
template <typename WImpl> using CompactWilsonClover = CompactWilsonCloverFermion<WImpl, CompactCloverHelpers<WImpl>>;
|
||||
template <typename WImpl> using CompactWilsonClover5D = CompactWilsonCloverFermion5D<WImpl, CompactCloverHelpers<WImpl>>;
|
||||
template <typename WImpl> using CompactWilsonExpClover = CompactWilsonCloverFermion<WImpl, CompactExpCloverHelpers<WImpl>>;
|
||||
|
||||
typedef CompactWilsonClover<WilsonImplD2> CompactWilsonCloverFermionD2;
|
||||
typedef CompactWilsonClover<WilsonImplF> CompactWilsonCloverFermionF;
|
||||
typedef CompactWilsonClover<WilsonImplD> CompactWilsonCloverFermionD;
|
||||
|
||||
typedef CompactWilsonClover5D<WilsonImplD2> CompactWilsonCloverFermion5DD2;
|
||||
typedef CompactWilsonClover5D<WilsonImplF> CompactWilsonCloverFermion5DF;
|
||||
typedef CompactWilsonClover5D<WilsonImplD> CompactWilsonCloverFermion5DD;
|
||||
|
||||
typedef CompactWilsonExpClover<WilsonImplD2> CompactWilsonExpCloverFermionD2;
|
||||
typedef CompactWilsonExpClover<WilsonImplF> CompactWilsonExpCloverFermionF;
|
||||
typedef CompactWilsonExpClover<WilsonImplD> CompactWilsonExpCloverFermionD;
|
||||
|
@ -91,13 +91,13 @@ public:
|
||||
virtual void Mdag (const FermionField &in, FermionField &out){assert(0);};
|
||||
|
||||
// half checkerboard operations; leave unimplemented as abstract for now
|
||||
virtual void Meooe (const FermionField &in, FermionField &out){assert(0);};
|
||||
virtual void Mooee (const FermionField &in, FermionField &out){assert(0);};
|
||||
virtual void MooeeInv (const FermionField &in, FermionField &out){assert(0);};
|
||||
virtual void Meooe (const FermionField &in, FermionField &out);
|
||||
virtual void Mooee (const FermionField &in, FermionField &out);
|
||||
virtual void MooeeInv (const FermionField &in, FermionField &out);
|
||||
|
||||
virtual void MeooeDag (const FermionField &in, FermionField &out){assert(0);};
|
||||
virtual void MooeeDag (const FermionField &in, FermionField &out){assert(0);};
|
||||
virtual void MooeeInvDag (const FermionField &in, FermionField &out){assert(0);};
|
||||
virtual void MeooeDag (const FermionField &in, FermionField &out);
|
||||
virtual void MooeeDag (const FermionField &in, FermionField &out);
|
||||
virtual void MooeeInvDag (const FermionField &in, FermionField &out);
|
||||
virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp){assert(0);}; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
|
||||
virtual void MdirAll(const FermionField &in, std::vector<FermionField> &out){assert(0);}; // case by case Wilson, Clover, Cayley, ContFrac, PartFrac
|
||||
|
||||
|
@ -0,0 +1,376 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/qcd/action/fermion/CompactWilsonCloverFermion5DImplementation.h
|
||||
|
||||
Copyright (C) 2017 - 2025
|
||||
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
||||
Author: Daniel Richtmann <daniel.richtmann@gmail.com>
|
||||
Author: Christoph Lehner <christoph@lhnr.de>
|
||||
|
||||
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 <Grid/Grid.h>
|
||||
#include <Grid/qcd/spin/Dirac.h>
|
||||
#include <Grid/qcd/action/fermion/CompactWilsonCloverFermion5D.h>
|
||||
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
template<class Impl, class CloverHelpers>
|
||||
CompactWilsonCloverFermion5D<Impl, CloverHelpers>::CompactWilsonCloverFermion5D(GaugeField& _Umu,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
const RealD _mass,
|
||||
const RealD _csw_r,
|
||||
const RealD _csw_t,
|
||||
const RealD _cF,
|
||||
const ImplParams& impl_p)
|
||||
: WilsonBase(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid, _mass, impl_p)
|
||||
, csw_r(_csw_r)
|
||||
, csw_t(_csw_t)
|
||||
, cF(_cF)
|
||||
, fixedBoundaries(impl_p.boundary_phases[Nd-1] == 0.0)
|
||||
, Diagonal(&FourDimGrid), Triangle(&FourDimGrid)
|
||||
, DiagonalEven(&FourDimRedBlackGrid), TriangleEven(&FourDimRedBlackGrid)
|
||||
, DiagonalOdd(&FourDimRedBlackGrid), TriangleOdd(&FourDimRedBlackGrid)
|
||||
, DiagonalInv(&FourDimGrid), TriangleInv(&FourDimGrid)
|
||||
, DiagonalInvEven(&FourDimRedBlackGrid), TriangleInvEven(&FourDimRedBlackGrid)
|
||||
, DiagonalInvOdd(&FourDimRedBlackGrid), TriangleInvOdd(&FourDimRedBlackGrid)
|
||||
, Tmp(&FiveDimGrid)
|
||||
, BoundaryMask(&FiveDimGrid)
|
||||
, BoundaryMaskEven(&FiveDimRedBlackGrid), BoundaryMaskOdd(&FiveDimRedBlackGrid)
|
||||
{
|
||||
assert(Nd == 4 && Nc == 3 && Ns == 4 && Impl::Dimension == 3);
|
||||
|
||||
csw_r *= 0.5;
|
||||
csw_t *= 0.5;
|
||||
//if (clover_anisotropy.isAnisotropic)
|
||||
// csw_r /= clover_anisotropy.xi_0;
|
||||
|
||||
ImportGauge(_Umu);
|
||||
if (fixedBoundaries) {
|
||||
this->BoundaryMaskEven.Checkerboard() = Even;
|
||||
this->BoundaryMaskOdd.Checkerboard() = Odd;
|
||||
CompactHelpers::SetupMasks(this->BoundaryMask, this->BoundaryMaskEven, this->BoundaryMaskOdd);
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::Dhop(const FermionField& in, FermionField& out, int dag) {
|
||||
WilsonBase::Dhop(in, out, dag);
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::DhopOE(const FermionField& in, FermionField& out, int dag) {
|
||||
WilsonBase::DhopOE(in, out, dag);
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::DhopEO(const FermionField& in, FermionField& out, int dag) {
|
||||
WilsonBase::DhopEO(in, out, dag);
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) {
|
||||
WilsonBase::DhopDir(in, out, dir, disp);
|
||||
if(this->fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::DhopDirAll(const FermionField& in, std::vector<FermionField>& out) {
|
||||
WilsonBase::DhopDirAll(in, out);
|
||||
if(this->fixedBoundaries) {
|
||||
for(auto& o : out) ApplyBoundaryMask(o);
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::M(const FermionField& in, FermionField& out) {
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
WilsonBase::Dhop(in, out, DaggerNo); // call base to save applying bc
|
||||
Mooee(in, Tmp);
|
||||
axpy(out, 1.0, out, Tmp);
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::Mdag(const FermionField& in, FermionField& out) {
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
WilsonBase::Dhop(in, out, DaggerYes); // call base to save applying bc
|
||||
MooeeDag(in, Tmp);
|
||||
axpy(out, 1.0, out, Tmp);
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::Meooe(const FermionField& in, FermionField& out) {
|
||||
WilsonBase::Meooe(in, out);
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MeooeDag(const FermionField& in, FermionField& out) {
|
||||
WilsonBase::MeooeDag(in, out);
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::Mooee(const FermionField& in, FermionField& out) {
|
||||
if(in.Grid()->_isCheckerBoarded) {
|
||||
if(in.Checkerboard() == Odd) {
|
||||
MooeeInternal(in, out, DiagonalOdd, TriangleOdd);
|
||||
} else {
|
||||
MooeeInternal(in, out, DiagonalEven, TriangleEven);
|
||||
}
|
||||
} else {
|
||||
MooeeInternal(in, out, Diagonal, Triangle);
|
||||
}
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MooeeDag(const FermionField& in, FermionField& out) {
|
||||
Mooee(in, out); // blocks are hermitian
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MooeeInv(const FermionField& in, FermionField& out) {
|
||||
if(in.Grid()->_isCheckerBoarded) {
|
||||
if(in.Checkerboard() == Odd) {
|
||||
MooeeInternal(in, out, DiagonalInvOdd, TriangleInvOdd);
|
||||
} else {
|
||||
MooeeInternal(in, out, DiagonalInvEven, TriangleInvEven);
|
||||
}
|
||||
} else {
|
||||
MooeeInternal(in, out, DiagonalInv, TriangleInv);
|
||||
}
|
||||
if(fixedBoundaries) ApplyBoundaryMask(out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MooeeInvDag(const FermionField& in, FermionField& out) {
|
||||
MooeeInv(in, out); // blocks are hermitian
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::Mdir(const FermionField& in, FermionField& out, int dir, int disp) {
|
||||
DhopDir(in, out, dir, disp);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MdirAll(const FermionField& in, std::vector<FermionField>& out) {
|
||||
DhopDirAll(in, out);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) {
|
||||
assert(!fixedBoundaries); // TODO check for changes required for open bc
|
||||
|
||||
// NOTE: code copied from original clover term
|
||||
conformable(X.Grid(), Y.Grid());
|
||||
conformable(X.Grid(), force.Grid());
|
||||
GaugeLinkField force_mu(force.Grid()), lambda(force.Grid());
|
||||
GaugeField clover_force(force.Grid());
|
||||
PropagatorField Lambda(force.Grid());
|
||||
|
||||
// Guido: Here we are hitting some performance issues:
|
||||
// need to extract the components of the DoubledGaugeField
|
||||
// for each call
|
||||
// Possible solution
|
||||
// Create a vector object to store them? (cons: wasting space)
|
||||
std::vector<GaugeLinkField> U(Nd, this->Umu.Grid());
|
||||
|
||||
Impl::extractLinkField(U, this->Umu);
|
||||
|
||||
force = Zero();
|
||||
// Derivative of the Wilson hopping term
|
||||
this->DhopDeriv(force, X, Y, dag);
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// Clover term derivative
|
||||
///////////////////////////////////////////////////////////
|
||||
Impl::outerProductImpl(Lambda, X, Y);
|
||||
//std::cout << "Lambda:" << Lambda << std::endl;
|
||||
|
||||
Gamma::Algebra sigma[] = {
|
||||
Gamma::Algebra::SigmaXY,
|
||||
Gamma::Algebra::SigmaXZ,
|
||||
Gamma::Algebra::SigmaXT,
|
||||
Gamma::Algebra::MinusSigmaXY,
|
||||
Gamma::Algebra::SigmaYZ,
|
||||
Gamma::Algebra::SigmaYT,
|
||||
Gamma::Algebra::MinusSigmaXZ,
|
||||
Gamma::Algebra::MinusSigmaYZ,
|
||||
Gamma::Algebra::SigmaZT,
|
||||
Gamma::Algebra::MinusSigmaXT,
|
||||
Gamma::Algebra::MinusSigmaYT,
|
||||
Gamma::Algebra::MinusSigmaZT};
|
||||
|
||||
/*
|
||||
sigma_{\mu \nu}=
|
||||
| 0 sigma[0] sigma[1] sigma[2] |
|
||||
| sigma[3] 0 sigma[4] sigma[5] |
|
||||
| sigma[6] sigma[7] 0 sigma[8] |
|
||||
| sigma[9] sigma[10] sigma[11] 0 |
|
||||
*/
|
||||
|
||||
int count = 0;
|
||||
clover_force = Zero();
|
||||
for (int mu = 0; mu < 4; mu++)
|
||||
{
|
||||
force_mu = Zero();
|
||||
for (int nu = 0; nu < 4; nu++)
|
||||
{
|
||||
if (mu == nu)
|
||||
continue;
|
||||
|
||||
RealD factor;
|
||||
if (nu == 4 || mu == 4)
|
||||
{
|
||||
factor = 2.0 * csw_t;
|
||||
}
|
||||
else
|
||||
{
|
||||
factor = 2.0 * csw_r;
|
||||
}
|
||||
PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked
|
||||
Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok
|
||||
force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu); // checked
|
||||
count++;
|
||||
}
|
||||
|
||||
pokeLorentz(clover_force, U[mu] * force_mu, mu);
|
||||
}
|
||||
//clover_force *= csw;
|
||||
force += clover_force;
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
|
||||
assert(0);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) {
|
||||
assert(0);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::MooeeInternal(const FermionField& in,
|
||||
FermionField& out,
|
||||
const CloverDiagonalField& diagonal,
|
||||
const CloverTriangleField& triangle) {
|
||||
assert(in.Checkerboard() == Odd || in.Checkerboard() == Even);
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
conformable(in, out);
|
||||
CompactHelpers::MooeeKernel(diagonal.oSites(), this->Ls, in, out, diagonal, triangle);
|
||||
}
|
||||
|
||||
template<class Impl, class CloverHelpers>
|
||||
void CompactWilsonCloverFermion5D<Impl, CloverHelpers>::ImportGauge(const GaugeField& _Umu) {
|
||||
// NOTE: parts copied from original implementation
|
||||
|
||||
// Import gauge into base class
|
||||
double t0 = usecond();
|
||||
WilsonBase::ImportGauge(_Umu); // NOTE: called here and in wilson constructor -> performed twice, but can't avoid that
|
||||
|
||||
// Initialize temporary variables
|
||||
double t1 = usecond();
|
||||
conformable(_Umu.Grid(), this->GaugeGrid());
|
||||
GridBase* grid = _Umu.Grid();
|
||||
typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid);
|
||||
CloverField TmpOriginal(grid);
|
||||
CloverField TmpInverse(grid);
|
||||
|
||||
// Compute the field strength terms mu>nu
|
||||
double t2 = usecond();
|
||||
WilsonLoops<Impl>::FieldStrength(Bx, _Umu, Zdir, Ydir);
|
||||
WilsonLoops<Impl>::FieldStrength(By, _Umu, Zdir, Xdir);
|
||||
WilsonLoops<Impl>::FieldStrength(Bz, _Umu, Ydir, Xdir);
|
||||
WilsonLoops<Impl>::FieldStrength(Ex, _Umu, Tdir, Xdir);
|
||||
WilsonLoops<Impl>::FieldStrength(Ey, _Umu, Tdir, Ydir);
|
||||
WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
|
||||
|
||||
// Compute the Clover Operator acting on Colour and Spin
|
||||
// multiply here by the clover coefficients for the anisotropy
|
||||
double t3 = usecond();
|
||||
TmpOriginal = Helpers::fillCloverYZ(Bx) * csw_r;
|
||||
TmpOriginal += Helpers::fillCloverXZ(By) * csw_r;
|
||||
TmpOriginal += Helpers::fillCloverXY(Bz) * csw_r;
|
||||
TmpOriginal += Helpers::fillCloverXT(Ex) * csw_t;
|
||||
TmpOriginal += Helpers::fillCloverYT(Ey) * csw_t;
|
||||
TmpOriginal += Helpers::fillCloverZT(Ez) * csw_t;
|
||||
|
||||
// Instantiate the clover term
|
||||
// - In case of the standard clover the mass term is added
|
||||
// - In case of the exponential clover the clover term is exponentiated
|
||||
double t4 = usecond();
|
||||
CloverHelpers::InstantiateClover(TmpOriginal, TmpInverse, csw_t, 4.0 + this->M5 /*this->diag_mass*/);
|
||||
|
||||
// Convert the data layout of the clover term
|
||||
double t5 = usecond();
|
||||
CompactHelpers::ConvertLayout(TmpOriginal, Diagonal, Triangle);
|
||||
|
||||
// Modify the clover term at the temporal boundaries in case of open boundary conditions
|
||||
double t6 = usecond();
|
||||
if(fixedBoundaries) CompactHelpers::ModifyBoundaries(Diagonal, Triangle, csw_t, cF, 4.0 + this->M5 /*this->diag_mass*/);
|
||||
|
||||
// Invert the Clover term
|
||||
// In case of the exponential clover with (anti-)periodic boundary conditions exp(-Clover) saved
|
||||
// in TmpInverse can be used. In all other cases the clover term has to be explictly inverted.
|
||||
// TODO: For now this inversion is explictly done on the CPU
|
||||
double t7 = usecond();
|
||||
CloverHelpers::InvertClover(TmpInverse, Diagonal, Triangle, DiagonalInv, TriangleInv, fixedBoundaries);
|
||||
|
||||
// Fill the remaining clover fields
|
||||
double t8 = usecond();
|
||||
pickCheckerboard(Even, DiagonalEven, Diagonal);
|
||||
pickCheckerboard(Even, TriangleEven, Triangle);
|
||||
pickCheckerboard(Odd, DiagonalOdd, Diagonal);
|
||||
pickCheckerboard(Odd, TriangleOdd, Triangle);
|
||||
pickCheckerboard(Even, DiagonalInvEven, DiagonalInv);
|
||||
pickCheckerboard(Even, TriangleInvEven, TriangleInv);
|
||||
pickCheckerboard(Odd, DiagonalInvOdd, DiagonalInv);
|
||||
pickCheckerboard(Odd, TriangleInvOdd, TriangleInv);
|
||||
|
||||
// Report timings
|
||||
double t9 = usecond();
|
||||
|
||||
std::cout << GridLogDebug << "CompactWilsonCloverFermion5D::ImportGauge timings:" << std::endl;
|
||||
std::cout << GridLogDebug << "WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "allocations = " << (t2 - t1) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "field strength = " << (t3 - t2) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "fill clover = " << (t4 - t3) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "instantiate clover = " << (t5 - t4) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "convert layout = " << (t6 - t5) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "modify boundaries = " << (t7 - t6) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "invert clover = " << (t8 - t7) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "pick cbs = " << (t9 - t8) / 1e6 << std::endl;
|
||||
std::cout << GridLogDebug << "total = " << (t9 - t0) / 1e6 << std::endl;
|
||||
}
|
||||
|
||||
NAMESPACE_END(Grid);
|
@ -14,6 +14,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
Author: Vera Guelpers <V.M.Guelpers@soton.ac.uk>
|
||||
Author: Christoph Lehner <christoph@lhnr.de>
|
||||
|
||||
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
|
||||
@ -484,6 +485,54 @@ void WilsonFermion5D<Impl>::DW(const FermionField &in, FermionField &out,int dag
|
||||
Dhop(in,out,dag); // -0.5 is included
|
||||
axpy(out,4.0-M5,in,out);
|
||||
}
|
||||
template <class Impl>
|
||||
void WilsonFermion5D<Impl>::Meooe(const FermionField &in, FermionField &out)
|
||||
{
|
||||
if (in.Checkerboard() == Odd) {
|
||||
DhopEO(in, out, DaggerNo);
|
||||
} else {
|
||||
DhopOE(in, out, DaggerNo);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion5D<Impl>::MeooeDag(const FermionField &in, FermionField &out)
|
||||
{
|
||||
if (in.Checkerboard() == Odd) {
|
||||
DhopEO(in, out, DaggerYes);
|
||||
} else {
|
||||
DhopOE(in, out, DaggerYes);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion5D<Impl>::Mooee(const FermionField &in, FermionField &out)
|
||||
{
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
typename FermionField::scalar_type scal(4.0 + M5);
|
||||
out = scal * in;
|
||||
}
|
||||
|
||||
template <class Impl>
|
||||
void WilsonFermion5D<Impl>::MooeeDag(const FermionField &in, FermionField &out)
|
||||
{
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
Mooee(in, out);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::MooeeInv(const FermionField &in, FermionField &out)
|
||||
{
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
out = (1.0/(4.0 + M5))*in;
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::MooeeInvDag(const FermionField &in, FermionField &out)
|
||||
{
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
MooeeInv(in,out);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt_5d(FermionField &out,const FermionField &in, RealD mass,std::vector<double> twist)
|
||||
|
@ -0,0 +1,45 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/ qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation5D.cc.master
|
||||
|
||||
Copyright (C) 2017 - 2025
|
||||
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Guido Cossu <guido.cossu@ed.ac.uk>
|
||||
Author: Daniel Richtmann <daniel.richtmann@gmail.com>
|
||||
Author: Mattia Bruno <mattia.bruno@cern.ch>
|
||||
Author: Christoph Lehner <christoph@lhnr.de>
|
||||
|
||||
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 <Grid/Grid.h>
|
||||
#include <Grid/qcd/spin/Dirac.h>
|
||||
#include <Grid/qcd/action/fermion/CompactWilsonCloverFermion5D.h>
|
||||
#include <Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermion5DImplementation.h>
|
||||
#include <Grid/qcd/action/fermion/CloverHelpers.h>
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
#include "impl.h"
|
||||
template class CompactWilsonCloverFermion5D<IMPLEMENTATION, CompactCloverHelpers<IMPLEMENTATION>>;
|
||||
template class CompactWilsonCloverFermion5D<IMPLEMENTATION, CompactExpCloverHelpers<IMPLEMENTATION>>;
|
||||
|
||||
NAMESPACE_END(Grid);
|
@ -0,0 +1 @@
|
||||
../CompactWilsonCloverFermion5DInstantiation.cc.master
|
@ -0,0 +1 @@
|
||||
../CompactWilsonCloverFermion5DInstantiation.cc.master
|
@ -62,7 +62,7 @@ do
|
||||
done
|
||||
done
|
||||
|
||||
CC_LIST="CompactWilsonCloverFermionInstantiation"
|
||||
CC_LIST="CompactWilsonCloverFermionInstantiation CompactWilsonCloverFermion5DInstantiation"
|
||||
|
||||
for impl in $COMPACT_WILSON_IMPL_LIST
|
||||
do
|
||||
|
@ -76,27 +76,27 @@ public:
|
||||
return action;
|
||||
};
|
||||
|
||||
virtual void deriv(const GaugeField &Umu,GaugeField & dSdU) {
|
||||
virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
|
||||
//extend Ta to include Lorentz indexes
|
||||
RealD factor_p = c_plaq/RealD(Nc)*0.5;
|
||||
RealD factor_r = c_rect/RealD(Nc)*0.5;
|
||||
|
||||
GridBase *grid = Umu.Grid();
|
||||
GridBase *grid = U.Grid();
|
||||
|
||||
std::vector<GaugeLinkField> U (Nd,grid);
|
||||
std::vector<GaugeLinkField> Umu (Nd,grid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
Umu[mu] = PeekIndex<LorentzIndex>(U,mu);
|
||||
}
|
||||
std::vector<GaugeLinkField> RectStaple(Nd,grid), Staple(Nd,grid);
|
||||
WilsonLoops<Gimpl>::StapleAndRectStapleAll(Staple, RectStaple, U, workspace);
|
||||
WilsonLoops<Gimpl>::StapleAndRectStapleAll(Staple, RectStaple, Umu, workspace);
|
||||
|
||||
GaugeLinkField dSdU_mu(grid);
|
||||
GaugeLinkField staple(grid);
|
||||
|
||||
for (int mu=0; mu < Nd; mu++){
|
||||
dSdU_mu = Ta(U[mu]*Staple[mu])*factor_p;
|
||||
dSdU_mu = dSdU_mu + Ta(U[mu]*RectStaple[mu])*factor_r;
|
||||
|
||||
dSdU_mu = Ta(Umu[mu]*Staple[mu])*factor_p;
|
||||
dSdU_mu = dSdU_mu + Ta(Umu[mu]*RectStaple[mu])*factor_r;
|
||||
|
||||
PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
|
||||
}
|
||||
|
||||
|
@ -73,20 +73,23 @@ public:
|
||||
// extend Ta to include Lorentz indexes
|
||||
|
||||
RealD factor = 0.5 * beta / RealD(Nc);
|
||||
GridBase *grid = U.Grid();
|
||||
|
||||
GaugeLinkField Umu(U.Grid());
|
||||
GaugeLinkField dSdU_mu(U.Grid());
|
||||
GaugeLinkField dSdU_mu(grid);
|
||||
std::vector<GaugeLinkField> Umu(Nd, grid);
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
Umu[mu] = PeekIndex<LorentzIndex>(U, mu);
|
||||
}
|
||||
|
||||
Umu = PeekIndex<LorentzIndex>(U, mu);
|
||||
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
// Staple in direction mu
|
||||
WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
|
||||
dSdU_mu = Ta(Umu * dSdU_mu) * factor;
|
||||
|
||||
WilsonLoops<Gimpl>::Staple(dSdU_mu, Umu, mu);
|
||||
dSdU_mu = Ta(Umu[mu] * dSdU_mu) * factor;
|
||||
|
||||
PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
RealD beta;
|
||||
};
|
||||
|
@ -207,11 +207,14 @@ std::vector<RealD> WilsonFlowBase<Gimpl>::flowMeasureEnergyDensityCloverleaf(con
|
||||
}
|
||||
|
||||
template <class Gimpl>
|
||||
void WilsonFlowBase<Gimpl>::setDefaultMeasurements(int topq_meas_interval){
|
||||
addMeasurement(1, [](int step, RealD t, const typename Gimpl::GaugeField &U){
|
||||
void WilsonFlowBase<Gimpl>::setDefaultMeasurements(int meas_interval){
|
||||
addMeasurement(meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
|
||||
std::cout << GridLogMessage << "[WilsonFlow] Energy density (plaq) : " << step << " " << t << " " << energyDensityPlaquette(t,U) << std::endl;
|
||||
});
|
||||
addMeasurement(topq_meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
|
||||
addMeasurement(meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
|
||||
std::cout << GridLogMessage << "[WilsonFlow] Energy density (cloverleaf) : " << step << " " << t << " " << energyDensityCloverleaf(t,U) << std::endl;
|
||||
});
|
||||
addMeasurement(meas_interval, [](int step, RealD t, const typename Gimpl::GaugeField &U){
|
||||
std::cout << GridLogMessage << "[WilsonFlow] Top. charge : " << step << " " << WilsonLoops<Gimpl>::TopologicalCharge(U) << std::endl;
|
||||
});
|
||||
}
|
||||
|
@ -292,19 +292,21 @@ public:
|
||||
//////////////////////////////////////////////////
|
||||
// the sum over all nu-oriented staples for nu != mu on each site
|
||||
//////////////////////////////////////////////////
|
||||
static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
|
||||
static void Staple(GaugeMat &staple, const GaugeLorentz &U, int mu) {
|
||||
|
||||
GridBase *grid = Umu.Grid();
|
||||
|
||||
std::vector<GaugeMat> U(Nd, grid);
|
||||
std::vector<GaugeMat> Umu(Nd, U.grid());
|
||||
for (int d = 0; d < Nd; d++) {
|
||||
U[d] = PeekIndex<LorentzIndex>(Umu, d);
|
||||
Umu[d] = PeekIndex<LorentzIndex>(U, d);
|
||||
}
|
||||
Staple(staple, U, mu);
|
||||
Staple(staple, Umu, mu);
|
||||
}
|
||||
|
||||
static void Staple(GaugeMat &staple, const std::vector<GaugeMat> &U, int mu) {
|
||||
staple = Zero();
|
||||
static void Staple(GaugeMat &staple, const std::vector<GaugeMat> &Umu, int mu) {
|
||||
|
||||
autoView(staple_v, staple, AcceleratorWrite);
|
||||
accelerator_for(i, staple.Grid()->oSites(), Simd::Nsimd(), {
|
||||
staple_v[i] = Zero();
|
||||
});
|
||||
|
||||
for (int nu = 0; nu < Nd; nu++) {
|
||||
|
||||
@ -318,12 +320,12 @@ public:
|
||||
// |
|
||||
// __|
|
||||
//
|
||||
|
||||
|
||||
staple += Gimpl::ShiftStaple(
|
||||
Gimpl::CovShiftForward(
|
||||
U[nu], nu,
|
||||
Umu[nu], nu,
|
||||
Gimpl::CovShiftBackward(
|
||||
U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))),
|
||||
Umu[mu], mu, Gimpl::CovShiftIdentityBackward(Umu[nu], nu))),
|
||||
mu);
|
||||
|
||||
// __
|
||||
@ -333,8 +335,8 @@ public:
|
||||
//
|
||||
|
||||
staple += Gimpl::ShiftStaple(
|
||||
Gimpl::CovShiftBackward(U[nu], nu,
|
||||
Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
|
||||
Gimpl::CovShiftBackward(Umu[nu], nu,
|
||||
Gimpl::CovShiftBackward(Umu[mu], mu, Umu[nu])), mu);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
@ -25,13 +25,20 @@ directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
#if Nc == 3
|
||||
#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
|
||||
#include <Grid/qcd/smearing/JacobianAction.h>
|
||||
#endif
|
||||
|
||||
using namespace Grid;
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
#if Nc != 3
|
||||
#warning FTHMC2p1f will not work for Nc != 3
|
||||
std::cout << "This program will currently only work for Nc == 3." << std::endl;
|
||||
#else
|
||||
std::cout << std::setprecision(12);
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
@ -220,7 +227,6 @@ int main(int argc, char **argv)
|
||||
TheHMC.Run(SmearingPolicy); // for smearing
|
||||
|
||||
Grid_finalize();
|
||||
#endif
|
||||
} // main
|
||||
|
||||
|
||||
|
||||
|
@ -24,14 +24,22 @@ See the full license in the file "LICENSE" in the top level distribution
|
||||
directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
#if Nc == 3
|
||||
#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
|
||||
#include <Grid/qcd/smearing/JacobianAction.h>
|
||||
#endif
|
||||
|
||||
using namespace Grid;
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
#if Nc != 3
|
||||
#warning FTHMC2p1f_3GeV will not work for Nc != 3
|
||||
std::cout << "This program will currently only work for Nc == 3." << std::endl;
|
||||
#else
|
||||
std::cout << std::setprecision(12);
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
@ -220,6 +228,7 @@ int main(int argc, char **argv)
|
||||
TheHMC.Run(SmearingPolicy); // for smearing
|
||||
|
||||
Grid_finalize();
|
||||
#endif
|
||||
} // main
|
||||
|
||||
|
||||
|
@ -25,13 +25,20 @@ directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
#if Nc == 3
|
||||
#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
|
||||
#include <Grid/qcd/smearing/JacobianAction.h>
|
||||
#endif
|
||||
|
||||
using namespace Grid;
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
#if Nc != 3
|
||||
#warning HMC2p1f_3GeV will not work for Nc != 3
|
||||
std::cout << "This program will currently only work for Nc == 3." << std::endl;
|
||||
#else
|
||||
std::cout << std::setprecision(12);
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
@ -220,6 +227,7 @@ int main(int argc, char **argv)
|
||||
TheHMC.Run(SmearingPolicy); // for smearing
|
||||
|
||||
Grid_finalize();
|
||||
#endif
|
||||
} // main
|
||||
|
||||
|
||||
|
@ -151,7 +151,7 @@ AC_ARG_ENABLE([tracing],
|
||||
case ${ac_TRACING} in
|
||||
nvtx)
|
||||
AC_DEFINE([GRID_TRACING_NVTX],[1],[use NVTX])
|
||||
LIBS="${LIBS} -lnvToolsExt64_1"
|
||||
LIBS="${LIBS} -lnvToolsExt"
|
||||
;;
|
||||
roctx)
|
||||
AC_DEFINE([GRID_TRACING_ROCTX],[1],[use ROCTX])
|
||||
|
@ -93,10 +93,13 @@ int main(int argc, char ** argv)
|
||||
Real coeff = (width*width) / Real(4*Iterations);
|
||||
|
||||
chi=kronecker;
|
||||
|
||||
// chi = (1-p^2/2N)^N kronecker
|
||||
for(int n = 0; n < Iterations; ++n) {
|
||||
Laplacian.M(chi,psi);
|
||||
chi = chi - coeff*psi;
|
||||
RealD n2 = norm2(chi);
|
||||
chi = chi * (1.0/std::sqrt(n2));
|
||||
}
|
||||
|
||||
std::cout << " Wuppertal smeared operator is chi = \n" << chi <<std::endl;
|
||||
|
14
tests/lanczos/LanParams.xml
Normal file
14
tests/lanczos/LanParams.xml
Normal file
@ -0,0 +1,14 @@
|
||||
<?xml version="1.0"?>
|
||||
<grid>
|
||||
<LanczosParameters>
|
||||
<mass>0.00107</mass>
|
||||
<M5>1.8</M5>
|
||||
<Ls>48</Ls>
|
||||
<Nstop>10</Nstop>
|
||||
<Nk>15</Nk>
|
||||
<Np>85</Np>
|
||||
<ChebyLow>0.003</ChebyLow>
|
||||
<ChebyHigh>60</ChebyHigh>
|
||||
<ChebyOrder>201</ChebyOrder>
|
||||
</LanczosParameters>
|
||||
</grid>
|
346
tests/lanczos/Test_dwf_G5R5.cc
Normal file
346
tests/lanczos/Test_dwf_G5R5.cc
Normal file
@ -0,0 +1,346 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_G5R5.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Chulwoo Jung <chulwoo@bnl.gov>
|
||||
From Duo and Bob's Chirality study
|
||||
|
||||
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 <Grid/Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
;
|
||||
|
||||
//typedef WilsonFermionD FermionOp;
|
||||
typedef DomainWallFermionD FermionOp;
|
||||
typedef typename DomainWallFermionD::FermionField FermionField;
|
||||
|
||||
|
||||
RealD AllZero(RealD x) { return 0.; }
|
||||
|
||||
namespace Grid {
|
||||
|
||||
struct LanczosParameters: Serializable {
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
|
||||
RealD, mass ,
|
||||
RealD, M5 ,
|
||||
Integer, Ls,
|
||||
Integer, Nstop,
|
||||
Integer, Nk,
|
||||
Integer, Np,
|
||||
RealD, ChebyLow,
|
||||
RealD, ChebyHigh,
|
||||
Integer, ChebyOrder)
|
||||
// Integer, StartTrajectory,
|
||||
// Integer, Trajectories, /* @brief Number of sweeps in this run */
|
||||
// bool, MetropolisTest,
|
||||
// Integer, NoMetropolisUntil,
|
||||
// std::string, StartingType,
|
||||
// Integer, SW,
|
||||
// RealD, Kappa,
|
||||
// IntegratorParameters, MD)
|
||||
|
||||
LanczosParameters() {
|
||||
////////////////////////////// Default values
|
||||
mass = 0;
|
||||
// MetropolisTest = true;
|
||||
// NoMetropolisUntil = 10;
|
||||
// StartTrajectory = 0;
|
||||
// SW = 2;
|
||||
// Trajectories = 10;
|
||||
// StartingType = "HotStart";
|
||||
/////////////////////////////////
|
||||
}
|
||||
|
||||
template <class ReaderClass >
|
||||
LanczosParameters(Reader<ReaderClass> & TheReader){
|
||||
initialize(TheReader);
|
||||
}
|
||||
|
||||
template < class ReaderClass >
|
||||
void initialize(Reader<ReaderClass> &TheReader){
|
||||
// std::cout << GridLogMessage << "Reading HMC\n";
|
||||
read(TheReader, "HMC", *this);
|
||||
}
|
||||
|
||||
|
||||
void print_parameters() const {
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Trajectories : " << Trajectories << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n";
|
||||
// MD.print_parameters();
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
int main(int argc, char** argv) {
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
LanczosParameters LanParams;
|
||||
#if 1
|
||||
{
|
||||
XmlReader HMCrd("LanParams.xml");
|
||||
read(HMCrd,"LanczosParameters",LanParams);
|
||||
}
|
||||
#else
|
||||
{
|
||||
LanParams.mass = mass;
|
||||
}
|
||||
#endif
|
||||
std::cout << GridLogMessage<< LanParams <<std::endl;
|
||||
{
|
||||
XmlWriter HMCwr("LanParams.xml.out");
|
||||
write(HMCwr,"LanczosParameters",LanParams);
|
||||
}
|
||||
|
||||
int Ls=16;
|
||||
RealD M5=1.8;
|
||||
RealD mass = -1.0;
|
||||
|
||||
mass=LanParams.mass;
|
||||
Ls=LanParams.Ls;
|
||||
M5=LanParams.M5;
|
||||
|
||||
GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
|
||||
GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridRedBlackCartesian* UrbGrid =
|
||||
SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
// GridCartesian* FGrid = UGrid;
|
||||
// GridRedBlackCartesian* FrbGrid = UrbGrid;
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid);
|
||||
// printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid, FrbGrid);
|
||||
|
||||
std::vector<int> seeds4({1, 2, 3, 4});
|
||||
std::vector<int> seeds5({5, 6, 7, 8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
FieldMetaData header;
|
||||
std::string file("./config");
|
||||
|
||||
int precision32 = 0;
|
||||
int tworow = 0;
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
|
||||
/*
|
||||
std::vector<LatticeColourMatrix> U(4, UGrid);
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
|
||||
}
|
||||
*/
|
||||
|
||||
int Nstop = 10;
|
||||
int Nk = 20;
|
||||
int Np = 80;
|
||||
Nstop=LanParams.Nstop;
|
||||
Nk=LanParams.Nk;
|
||||
Np=LanParams.Np;
|
||||
|
||||
int Nm = Nk + Np;
|
||||
int MaxIt = 10000;
|
||||
RealD resid = 1.0e-5;
|
||||
|
||||
|
||||
//while ( mass > - 5.0){
|
||||
FermionOp Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
MdagMLinearOperator<FermionOp,FermionField> HermOp(Ddwf); /// <-----
|
||||
// Gamma5HermitianLinearOperator <FermionOp,LatticeFermion> HermOp2(WilsonOperator); /// <-----
|
||||
Gamma5R5HermitianLinearOperator<FermionOp, LatticeFermion> G5R5Herm(Ddwf);
|
||||
// Gamma5R5HermitianLinearOperator
|
||||
std::vector<double> Coeffs{0, 1.};
|
||||
Polynomial<FermionField> PolyX(Coeffs);
|
||||
|
||||
Chebyshev<FermionField> Cheby(LanParams.ChebyLow,LanParams.ChebyHigh,LanParams.ChebyOrder);
|
||||
|
||||
FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
|
||||
PlainHermOp<FermionField> Op (HermOp);
|
||||
PlainHermOp<FermionField> Op2 (G5R5Herm);
|
||||
|
||||
ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby, Op, Nstop, Nk, Nm, resid, MaxIt);
|
||||
|
||||
std::vector<RealD> eval(Nm);
|
||||
FermionField src(FGrid);
|
||||
gaussian(RNG5, src);
|
||||
std::vector<FermionField> evec(Nm, FGrid);
|
||||
for (int i = 0; i < 1; i++) {
|
||||
std::cout << i << " / " << Nm << " grid pointer " << evec[i].Grid()
|
||||
<< std::endl;
|
||||
};
|
||||
|
||||
int Nconv;
|
||||
IRL.calc(eval, evec, src, Nconv);
|
||||
|
||||
std::cout << mass <<" : " << eval << std::endl;
|
||||
|
||||
#if 0
|
||||
Gamma g5(Gamma::Algebra::Gamma5) ;
|
||||
ComplexD dot;
|
||||
FermionField tmp(FGrid);
|
||||
// RealD eMe,eMMe;
|
||||
for (int i = 0; i < Nstop ; i++) {
|
||||
// tmp = g5*evec[i];
|
||||
dot = innerProduct(evec[i],evec[i]);
|
||||
// G5R5(tmp,evec[i]);
|
||||
G5R5Herm.HermOpAndNorm(evec[i],tmp,eMe,eMMe);
|
||||
std::cout <<"Norm "<<M5<<" "<< mass << " : " << i << " " << real(dot) << " " << imag(dot) << " "<< eMe << " " <<eMMe<< std::endl ;
|
||||
for (int j = 0; j < Nstop ; j++) {
|
||||
dot = innerProduct(tmp,evec[j]);
|
||||
std::cout <<"G5R5 "<<M5<<" "<< mass << " : " << i << " " <<j<<" " << real(dot) << " " << imag(dot) << std::endl ;
|
||||
}
|
||||
}
|
||||
// src = evec[0]+evec[1]+evec[2];
|
||||
// mass += -0.1;
|
||||
#endif
|
||||
|
||||
//**********************************************************************
|
||||
//orthogonalization
|
||||
//calculat the matrix
|
||||
cout << "Start orthogonalization " << endl;
|
||||
cout << "calculate the matrix element" << endl;
|
||||
vector<LatticeFermion> G5R5Mevec(Nconv, FGrid);
|
||||
vector<LatticeFermion> finalevec(Nconv, FGrid);
|
||||
vector<RealD> eMe(Nconv), eMMe(Nconv);
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
G5R5Herm.HermOpAndNorm(evec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
|
||||
}
|
||||
cout << "Re<evec, G5R5M(evec)>: " << endl;
|
||||
cout << eMe << endl;
|
||||
cout << "<G5R5M(evec), G5R5M(evec)>" << endl;
|
||||
cout << eMMe << endl;
|
||||
vector<vector<ComplexD>> VevecG5R5Mevec(Nconv);
|
||||
Eigen::MatrixXcd evecG5R5Mevec = Eigen::MatrixXcd::Zero(Nconv, Nconv);
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
VevecG5R5Mevec[i].resize(Nconv);
|
||||
for(int j = 0; j < Nconv; j++){
|
||||
VevecG5R5Mevec[i][j] = innerProduct(evec[i], G5R5Mevec[j]);
|
||||
evecG5R5Mevec(i, j) = VevecG5R5Mevec[i][j];
|
||||
}
|
||||
}
|
||||
//calculate eigenvector
|
||||
cout << "Eigen solver" << endl;
|
||||
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver(evecG5R5Mevec);
|
||||
vector<RealD> eigeneval(Nconv);
|
||||
vector<vector<ComplexD>> eigenevec(Nconv);
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
eigeneval[i] = eigensolver.eigenvalues()[i];
|
||||
eigenevec[i].resize(Nconv);
|
||||
for(int j = 0; j < Nconv; j++){
|
||||
eigenevec[i][j] = eigensolver.eigenvectors()(i, j);
|
||||
}
|
||||
}
|
||||
//rotation
|
||||
cout << "Do rotation" << endl;
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
finalevec[i] = finalevec[i] - finalevec[i];
|
||||
for(int j = 0; j < Nconv; j++){
|
||||
finalevec[i] = eigenevec[j][i]*evec[j] + finalevec[i];
|
||||
}
|
||||
}
|
||||
//normalize again;
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
RealD tmp_RealD = norm2(finalevec[i]);
|
||||
tmp_RealD = 1./pow(tmp_RealD, 0.5);
|
||||
finalevec[i] = finalevec[i]*tmp_RealD;
|
||||
}
|
||||
|
||||
//check
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
G5R5Herm.HermOpAndNorm(finalevec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
|
||||
}
|
||||
|
||||
//**********************************************************************
|
||||
//sort the eigenvectors
|
||||
vector<LatticeFermion> finalevec_copy(Nconv, FGrid);
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
finalevec_copy[i] = finalevec[i];
|
||||
}
|
||||
vector<RealD> eMe_copy(eMe);
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
eMe[i] = fabs(eMe[i]);
|
||||
eMe_copy[i] = eMe[i];
|
||||
}
|
||||
sort(eMe_copy.begin(), eMe_copy.end());
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
for(int j = 0; j < Nconv; j++){
|
||||
if(eMe[j] == eMe_copy[i]){
|
||||
finalevec[i] = finalevec_copy[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
G5R5Herm.HermOpAndNorm(finalevec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
|
||||
}
|
||||
cout << "Re<evec, G5R5M(evec)>: " << endl;
|
||||
cout << eMe << endl;
|
||||
cout << "<G5R5M(evec), G5R5M(evec)>" << endl;
|
||||
cout << eMMe << endl;
|
||||
|
||||
|
||||
// vector<LatticeFermion> finalevec(Nconv, FGrid);
|
||||
// temporary, until doing rotation
|
||||
// for(int i = 0; i < Nconv; i++)
|
||||
// finalevec[i]=evec[i];
|
||||
//**********************************************************************
|
||||
//calculate chirality matrix
|
||||
vector<LatticeFermion> G5evec(Nconv, FGrid);
|
||||
vector<vector<ComplexD>> chiral_matrix(Nconv);
|
||||
vector<vector<RealD>> chiral_matrix_real(Nconv);
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
// G5evec[i] = G5evec[i] - G5evec[i];
|
||||
G5evec[i] = Zero();
|
||||
for(int j = 0; j < Ls/2; j++){
|
||||
axpby_ssp(G5evec[i], 1., finalevec[i], 0., G5evec[i], j, j);
|
||||
}
|
||||
for(int j = Ls/2; j < Ls; j++){
|
||||
axpby_ssp(G5evec[i], -1., finalevec[i], 0., G5evec[i], j, j);
|
||||
}
|
||||
}
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
chiral_matrix_real[i].resize(Nconv);
|
||||
chiral_matrix[i].resize(Nconv);
|
||||
for(int j = 0; j < Nconv; j++){
|
||||
chiral_matrix[i][j] = innerProduct(finalevec[i], G5evec[j]);
|
||||
chiral_matrix_real[i][j] = abs(chiral_matrix[i][j]);
|
||||
std::cout <<" chiral_matrix_real "<<i<<" "<<j<<" "<< chiral_matrix_real[i][j] << std::endl;
|
||||
}
|
||||
}
|
||||
for(int i = 0; i < Nconv; i++){
|
||||
if(chiral_matrix[i][i].real() < 0.){
|
||||
chiral_matrix_real[i][i] = -1. * chiral_matrix_real[i][i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
278
tests/lanczos/Test_wilson_DWFKernel.cc
Normal file
278
tests/lanczos/Test_wilson_DWFKernel.cc
Normal file
@ -0,0 +1,278 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_lanczos.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Chulwoo Jung <chulwoo@bnl.gov>
|
||||
|
||||
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 <Grid/Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
;
|
||||
|
||||
typedef WilsonFermionD FermionOp;
|
||||
typedef typename WilsonFermionD::FermionField FermionField;
|
||||
|
||||
|
||||
RealD AllZero(RealD x) { return 0.; }
|
||||
|
||||
namespace Grid {
|
||||
|
||||
#if 0
|
||||
template<typename Field>
|
||||
class RationalHermOp : public LinearFunction<Field> {
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
// OperatorFunction<Field> & _poly;
|
||||
LinearOperatorBase<Field> &_Linop;
|
||||
RealD _massDen, _massNum;
|
||||
|
||||
FunctionHermOp(LinearOperatorBase<Field>& linop, RealD massDen,RealD massNum)
|
||||
: _Linop(linop) ,_massDen(massDen),_massNum(massNum) {};
|
||||
|
||||
void operator()(const Field& in, Field& out) {
|
||||
// _poly(_Linop,in,out);
|
||||
}
|
||||
};
|
||||
#endif
|
||||
|
||||
template<class Matrix,class Field>
|
||||
class InvG5LinearOperator : public LinearOperatorBase<Field> {
|
||||
Matrix &_Mat;
|
||||
RealD _num;
|
||||
RealD _Tol;
|
||||
Integer _MaxIt;
|
||||
Gamma g5;
|
||||
|
||||
public:
|
||||
InvG5LinearOperator(Matrix &Mat,RealD num): _Mat(Mat),_num(num), _Tol(1e-12),_MaxIt(10000), g5(Gamma::Algebra::Gamma5) {};
|
||||
|
||||
// Support for coarsening to a multigrid
|
||||
void OpDiag (const Field &in, Field &out) {
|
||||
assert(0);
|
||||
_Mat.Mdiag(in,out);
|
||||
}
|
||||
void OpDir (const Field &in, Field &out,int dir,int disp) {
|
||||
assert(0);
|
||||
_Mat.Mdir(in,out,dir,disp);
|
||||
}
|
||||
void OpDirAll (const Field &in, std::vector<Field> &out){
|
||||
assert(0);
|
||||
_Mat.MdirAll(in,out);
|
||||
};
|
||||
void Op (const Field &in, Field &out){
|
||||
assert(0);
|
||||
_Mat.M(in,out);
|
||||
}
|
||||
void AdjOp (const Field &in, Field &out){
|
||||
assert(0);
|
||||
_Mat.Mdag(in,out);
|
||||
}
|
||||
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
|
||||
HermOp(in,out);
|
||||
ComplexD dot = innerProduct(in,out);
|
||||
n1=real(dot);
|
||||
n2=norm2(out);
|
||||
}
|
||||
void HermOp(const Field &in, Field &out){
|
||||
Field tmp(in.Grid());
|
||||
MdagMLinearOperator<Matrix,Field> denom(_Mat);
|
||||
ConjugateGradient<Field> CG(_Tol,_MaxIt);
|
||||
_Mat.M(in,tmp);
|
||||
tmp += _num*in;
|
||||
_Mat.Mdag(tmp,out);
|
||||
CG(denom,out,tmp);
|
||||
out = g5*tmp;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
struct LanczosParameters: Serializable {
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
|
||||
RealD, mass ,
|
||||
RealD, resid,
|
||||
RealD, ChebyLow,
|
||||
RealD, ChebyHigh,
|
||||
Integer, ChebyOrder)
|
||||
// Integer, StartTrajectory,
|
||||
// Integer, Trajectories, /* @brief Number of sweeps in this run */
|
||||
// bool, MetropolisTest,
|
||||
// Integer, NoMetropolisUntil,
|
||||
// std::string, StartingType,
|
||||
// Integer, SW,
|
||||
// RealD, Kappa,
|
||||
// IntegratorParameters, MD)
|
||||
|
||||
LanczosParameters() {
|
||||
////////////////////////////// Default values
|
||||
mass = 0;
|
||||
// MetropolisTest = true;
|
||||
// NoMetropolisUntil = 10;
|
||||
// StartTrajectory = 0;
|
||||
// SW = 2;
|
||||
// Trajectories = 10;
|
||||
// StartingType = "HotStart";
|
||||
/////////////////////////////////
|
||||
}
|
||||
|
||||
template <class ReaderClass >
|
||||
LanczosParameters(Reader<ReaderClass> & TheReader){
|
||||
initialize(TheReader);
|
||||
}
|
||||
|
||||
template < class ReaderClass >
|
||||
void initialize(Reader<ReaderClass> &TheReader){
|
||||
// std::cout << GridLogMessage << "Reading HMC\n";
|
||||
read(TheReader, "HMC", *this);
|
||||
}
|
||||
|
||||
|
||||
void print_parameters() const {
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Trajectories : " << Trajectories << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n";
|
||||
// MD.print_parameters();
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
int main(int argc, char** argv) {
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
|
||||
GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridRedBlackCartesian* UrbGrid =
|
||||
SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridCartesian* FGrid = UGrid;
|
||||
GridRedBlackCartesian* FrbGrid = UrbGrid;
|
||||
// printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid, FrbGrid);
|
||||
|
||||
std::vector<int> seeds4({1, 2, 3, 4});
|
||||
std::vector<int> seeds5({5, 6, 7, 8});
|
||||
GridParallelRNG RNG5(FGrid);
|
||||
RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid);
|
||||
RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG RNG5rb(FrbGrid);
|
||||
RNG5.SeedFixedIntegers(seeds5);
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
// SU<Nc>::HotConfiguration(RNG4, Umu);
|
||||
|
||||
FieldMetaData header;
|
||||
std::string file("./config");
|
||||
|
||||
int precision32 = 0;
|
||||
int tworow = 0;
|
||||
// NerscIO::writeConfiguration(Umu,file,tworow,precision32);
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
|
||||
/*
|
||||
std::vector<LatticeColourMatrix> U(4, UGrid);
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
|
||||
}
|
||||
*/
|
||||
|
||||
int Nstop = 5;
|
||||
int Nk = 10;
|
||||
int Np = 90;
|
||||
int Nm = Nk + Np;
|
||||
int MaxIt = 10000;
|
||||
RealD resid = 1.0e-5;
|
||||
|
||||
RealD mass = -1.0;
|
||||
|
||||
LanczosParameters LanParams;
|
||||
#if 1
|
||||
{
|
||||
XmlReader HMCrd("LanParams.xml");
|
||||
read(HMCrd,"LanczosParameters",LanParams);
|
||||
}
|
||||
#else
|
||||
{
|
||||
LanParams.mass = mass;
|
||||
}
|
||||
#endif
|
||||
std::cout << GridLogMessage<< LanParams <<std::endl;
|
||||
{
|
||||
XmlWriter HMCwr("LanParams.xml.out");
|
||||
write(HMCwr,"LanczosParameters",LanParams);
|
||||
}
|
||||
|
||||
mass=LanParams.mass;
|
||||
resid=LanParams.resid;
|
||||
|
||||
|
||||
while ( mass > - 5.0){
|
||||
FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,2.+mass);
|
||||
InvG5LinearOperator<FermionOp,LatticeFermion> HermOp(WilsonOperator,-2.); /// <-----
|
||||
//SchurDiagTwoOperator<FermionOp,FermionField> HermOp(WilsonOperator);
|
||||
// Gamma5HermitianLinearOperator <FermionOp,LatticeFermion> HermOp2(WilsonOperator); /// <-----
|
||||
|
||||
std::vector<double> Coeffs{0, 0, 1.};
|
||||
Polynomial<FermionField> PolyX(Coeffs);
|
||||
Chebyshev<FermionField> Cheby(LanParams.ChebyLow,LanParams.ChebyHigh,LanParams.ChebyOrder);
|
||||
|
||||
FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
|
||||
// InvHermOp<FermionField> Op(WilsonOperator,HermOp);
|
||||
PlainHermOp<FermionField> Op (HermOp);
|
||||
// PlainHermOp<FermionField> Op2 (HermOp2);
|
||||
|
||||
ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby, Op, Nstop, Nk, Nm, resid, MaxIt);
|
||||
|
||||
std::vector<RealD> eval(Nm);
|
||||
FermionField src(FGrid);
|
||||
gaussian(RNG5, src);
|
||||
std::vector<FermionField> evec(Nm, FGrid);
|
||||
for (int i = 0; i < 1; i++) {
|
||||
std::cout << i << " / " << Nm << " grid pointer " << evec[i].Grid()
|
||||
<< std::endl;
|
||||
};
|
||||
|
||||
int Nconv;
|
||||
IRL.calc(eval, evec, src, Nconv);
|
||||
|
||||
std::cout << mass <<" : " << eval << std::endl;
|
||||
|
||||
Gamma g5(Gamma::Algebra::Gamma5) ;
|
||||
ComplexD dot;
|
||||
FermionField tmp(FGrid);
|
||||
for (int i = 0; i < Nstop ; i++) {
|
||||
tmp = g5*evec[i];
|
||||
dot = innerProduct(tmp,evec[i]);
|
||||
std::cout << mass << " : " << eval[i] << " " << real(dot) << " " << imag(dot) << std::endl ;
|
||||
}
|
||||
src = evec[0]+evec[1]+evec[2];
|
||||
mass += -0.1;
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
211
tests/lanczos/Test_wilson_specflow.cc
Normal file
211
tests/lanczos/Test_wilson_specflow.cc
Normal file
@ -0,0 +1,211 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_dwf_lanczos.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Chulwoo Jung <chulwoo@bnl.gov>
|
||||
|
||||
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 <Grid/Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
;
|
||||
|
||||
typedef WilsonFermionD FermionOp;
|
||||
typedef typename WilsonFermionD::FermionField FermionField;
|
||||
|
||||
|
||||
RealD AllZero(RealD x) { return 0.; }
|
||||
|
||||
namespace Grid {
|
||||
|
||||
struct LanczosParameters: Serializable {
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(LanczosParameters,
|
||||
RealD, mass ,
|
||||
RealD, ChebyLow,
|
||||
RealD, ChebyHigh,
|
||||
Integer, ChebyOrder)
|
||||
// Integer, StartTrajectory,
|
||||
// Integer, Trajectories, /* @brief Number of sweeps in this run */
|
||||
// bool, MetropolisTest,
|
||||
// Integer, NoMetropolisUntil,
|
||||
// std::string, StartingType,
|
||||
// Integer, SW,
|
||||
// RealD, Kappa,
|
||||
// IntegratorParameters, MD)
|
||||
|
||||
LanczosParameters() {
|
||||
////////////////////////////// Default values
|
||||
mass = 0;
|
||||
// MetropolisTest = true;
|
||||
// NoMetropolisUntil = 10;
|
||||
// StartTrajectory = 0;
|
||||
// SW = 2;
|
||||
// Trajectories = 10;
|
||||
// StartingType = "HotStart";
|
||||
/////////////////////////////////
|
||||
}
|
||||
|
||||
template <class ReaderClass >
|
||||
LanczosParameters(Reader<ReaderClass> & TheReader){
|
||||
initialize(TheReader);
|
||||
}
|
||||
|
||||
template < class ReaderClass >
|
||||
void initialize(Reader<ReaderClass> &TheReader){
|
||||
// std::cout << GridLogMessage << "Reading HMC\n";
|
||||
read(TheReader, "HMC", *this);
|
||||
}
|
||||
|
||||
|
||||
void print_parameters() const {
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Trajectories : " << Trajectories << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Start trajectory : " << StartTrajectory << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Metropolis test (on/off): " << std::boolalpha << MetropolisTest << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Thermalization trajs : " << NoMetropolisUntil << "\n";
|
||||
// std::cout << GridLogMessage << "[HMC parameters] Starting type : " << StartingType << "\n";
|
||||
// MD.print_parameters();
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
int main(int argc, char** argv) {
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
|
||||
GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridRedBlackCartesian* UrbGrid =
|
||||
SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridCartesian* FGrid = UGrid;
|
||||
GridRedBlackCartesian* FrbGrid = UrbGrid;
|
||||
// printf("UGrid=%p UrbGrid=%p FGrid=%p FrbGrid=%p\n", UGrid, UrbGrid, FGrid, FrbGrid);
|
||||
|
||||
std::vector<int> seeds4({1, 2, 3, 4});
|
||||
std::vector<int> seeds5({5, 6, 7, 8});
|
||||
GridParallelRNG RNG5(FGrid);
|
||||
RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid);
|
||||
RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG RNG5rb(FrbGrid);
|
||||
RNG5.SeedFixedIntegers(seeds5);
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
// SU<Nc>::HotConfiguration(RNG4, Umu);
|
||||
|
||||
FieldMetaData header;
|
||||
std::string file("./config");
|
||||
|
||||
int precision32 = 0;
|
||||
int tworow = 0;
|
||||
// NerscIO::writeConfiguration(Umu,file,tworow,precision32);
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
|
||||
/*
|
||||
std::vector<LatticeColourMatrix> U(4, UGrid);
|
||||
for (int mu = 0; mu < Nd; mu++) {
|
||||
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
|
||||
}
|
||||
*/
|
||||
|
||||
int Nstop = 10;
|
||||
int Nk = 20;
|
||||
int Np = 80;
|
||||
int Nm = Nk + Np;
|
||||
int MaxIt = 10000;
|
||||
RealD resid = 1.0e-5;
|
||||
|
||||
RealD mass = -1.0;
|
||||
|
||||
LanczosParameters LanParams;
|
||||
#if 1
|
||||
{
|
||||
XmlReader HMCrd("LanParams.xml");
|
||||
read(HMCrd,"LanczosParameters",LanParams);
|
||||
}
|
||||
#else
|
||||
{
|
||||
LanParams.mass = mass;
|
||||
}
|
||||
#endif
|
||||
std::cout << GridLogMessage<< LanParams <<std::endl;
|
||||
{
|
||||
XmlWriter HMCwr("LanParams.xml.out");
|
||||
write(HMCwr,"LanczosParameters",LanParams);
|
||||
}
|
||||
|
||||
mass=LanParams.mass;
|
||||
|
||||
|
||||
while ( mass > - 5.0){
|
||||
FermionOp WilsonOperator(Umu,*FGrid,*FrbGrid,mass);
|
||||
MdagMLinearOperator<FermionOp,FermionField> HermOp(WilsonOperator); /// <-----
|
||||
//SchurDiagTwoOperator<FermionOp,FermionField> HermOp(WilsonOperator);
|
||||
Gamma5HermitianLinearOperator <FermionOp,LatticeFermion> HermOp2(WilsonOperator); /// <-----
|
||||
|
||||
std::vector<double> Coeffs{0, 1.};
|
||||
Polynomial<FermionField> PolyX(Coeffs);
|
||||
// Chebyshev<FermionField> Cheby(0.5, 60., 31);
|
||||
// RealD, ChebyLow,
|
||||
// RealD, ChebyHigh,
|
||||
// Integer, ChebyOrder)
|
||||
|
||||
Chebyshev<FermionField> Cheby(LanParams.ChebyLow,LanParams.ChebyHigh,LanParams.ChebyOrder);
|
||||
|
||||
FunctionHermOp<FermionField> OpCheby(Cheby,HermOp);
|
||||
PlainHermOp<FermionField> Op (HermOp);
|
||||
PlainHermOp<FermionField> Op2 (HermOp2);
|
||||
|
||||
ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby, Op2, Nstop, Nk, Nm, resid, MaxIt);
|
||||
|
||||
std::vector<RealD> eval(Nm);
|
||||
FermionField src(FGrid);
|
||||
gaussian(RNG5, src);
|
||||
std::vector<FermionField> evec(Nm, FGrid);
|
||||
for (int i = 0; i < 1; i++) {
|
||||
std::cout << i << " / " << Nm << " grid pointer " << evec[i].Grid()
|
||||
<< std::endl;
|
||||
};
|
||||
|
||||
int Nconv;
|
||||
IRL.calc(eval, evec, src, Nconv);
|
||||
|
||||
std::cout << mass <<" : " << eval << std::endl;
|
||||
|
||||
Gamma g5(Gamma::Algebra::Gamma5) ;
|
||||
ComplexD dot;
|
||||
FermionField tmp(FGrid);
|
||||
for (int i = 0; i < Nstop ; i++) {
|
||||
tmp = g5*evec[i];
|
||||
dot = innerProduct(tmp,evec[i]);
|
||||
std::cout << mass << " : " << eval[i] << " " << real(dot) << " " << imag(dot) << std::endl ;
|
||||
}
|
||||
src = evec[0]+evec[1]+evec[2];
|
||||
mass += -0.1;
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
@ -33,8 +33,7 @@ namespace Grid{
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(WFParameters,
|
||||
int, steps,
|
||||
double, step_size,
|
||||
int, meas_interval,
|
||||
double, maxTau); // for the adaptive algorithm
|
||||
int, meas_interval);
|
||||
|
||||
|
||||
template <class ReaderClass >
|
||||
@ -86,7 +85,7 @@ int main(int argc, char **argv) {
|
||||
WFParameters WFPar(Reader);
|
||||
ConfParameters CPar(Reader);
|
||||
CheckpointerParameters CPPar(CPar.conf_prefix, CPar.rng_prefix);
|
||||
BinaryHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
|
||||
NerscHmcCheckpointer<PeriodicGimplR> CPBin(CPPar);
|
||||
|
||||
for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
|
||||
|
||||
@ -96,19 +95,13 @@ int main(int argc, char **argv) {
|
||||
std::cout << GridLogMessage << "Initial plaquette: "
|
||||
<< WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
|
||||
|
||||
int t=WFPar.maxTau;
|
||||
WilsonFlowAdaptive<PeriodicGimplR> WF(WFPar.step_size, WFPar.maxTau,
|
||||
1.0e-4,
|
||||
WilsonFlow<PeriodicGimplR> WF(WFPar.step_size, WFPar.steps,
|
||||
WFPar.meas_interval);
|
||||
|
||||
WF.smear(Uflow, Umu);
|
||||
|
||||
RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
|
||||
RealD WFlow_TC = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);
|
||||
RealD WFlow_T0 = WF.energyDensityPlaquette(t,Uflow);
|
||||
std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl;
|
||||
std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl;
|
||||
std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << std::endl;
|
||||
|
||||
std::cout<< GridLogMessage << " Admissibility check:\n";
|
||||
const double sp_adm = 0.067; // admissible threshold
|
||||
|
37
visualisation/CMakeLists.txt
Normal file
37
visualisation/CMakeLists.txt
Normal file
@ -0,0 +1,37 @@
|
||||
cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
|
||||
|
||||
project(GridViewer)
|
||||
|
||||
list(APPEND CMAKE_PREFIX_PATH "/Users/peterboyle/QCD/vtk/VTK-9.4.2-install/")
|
||||
|
||||
find_package(VTK COMPONENTS
|
||||
CommonColor
|
||||
CommonCore
|
||||
FiltersCore
|
||||
FiltersModeling
|
||||
IOImage
|
||||
IOFFMPEG
|
||||
InteractionStyle
|
||||
InteractionWidgets
|
||||
RenderingContextOpenGL2
|
||||
RenderingCore
|
||||
RenderingFreeType
|
||||
RenderingGL2PSOpenGL2
|
||||
RenderingOpenGL2
|
||||
)
|
||||
|
||||
if (NOT VTK_FOUND)
|
||||
message(FATAL_ERROR "GridViewer: Unable to find the VTK build folder.")
|
||||
endif()
|
||||
|
||||
# Prevent a "command line is too long" failure in Windows.
|
||||
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
|
||||
|
||||
add_executable(FieldDensityAnimate MACOSX_BUNDLE FieldDensityAnimate.cxx )
|
||||
target_link_libraries(FieldDensityAnimate PRIVATE ${VTK_LIBRARIES}
|
||||
)
|
||||
# vtk_module_autoinit is needed
|
||||
vtk_module_autoinit(
|
||||
TARGETS FieldDensityAnimate
|
||||
MODULES ${VTK_LIBRARIES}
|
||||
)
|
BIN
visualisation/E8_vs_Topo8.avi
Normal file
BIN
visualisation/E8_vs_Topo8.avi
Normal file
Binary file not shown.
285
visualisation/FieldDensity.py
Normal file
285
visualisation/FieldDensity.py
Normal file
@ -0,0 +1,285 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
# noinspection PyUnresolvedReferences
|
||||
import math
|
||||
import vtk
|
||||
import vtkmodules.vtkInteractionStyle
|
||||
# noinspection PyUnresolvedReferences
|
||||
import vtkmodules.vtkRenderingOpenGL2
|
||||
from vtkmodules.vtkCommonColor import vtkNamedColors
|
||||
from vtkmodules.vtkCommonCore import (
|
||||
VTK_VERSION_NUMBER,
|
||||
vtkVersion
|
||||
)
|
||||
from vtkmodules.vtkCommonCore import VTK_DOUBLE
|
||||
from vtkmodules.vtkCommonDataModel import vtkImageData
|
||||
from vtkmodules.vtkFiltersCore import (
|
||||
vtkMarchingCubes,
|
||||
vtkStripper
|
||||
)
|
||||
from vtkmodules.vtkFiltersModeling import vtkOutlineFilter
|
||||
from vtkmodules.vtkIOImage import (
|
||||
vtkMetaImageReader,
|
||||
vtkJPEGWriter,
|
||||
vtkPNGWriter
|
||||
)
|
||||
from vtkmodules.vtkRenderingCore import (
|
||||
vtkActor,
|
||||
vtkCamera,
|
||||
vtkPolyDataMapper,
|
||||
vtkProperty,
|
||||
vtkRenderWindow,
|
||||
vtkRenderWindowInteractor,
|
||||
vtkRenderer,
|
||||
vtkWindowToImageFilter
|
||||
)
|
||||
|
||||
|
||||
class vtkTimerCallback():
|
||||
def __init__(self, steps, imageData, iren):
|
||||
self.timer_count = 0
|
||||
self.steps = steps
|
||||
self.imageData = imageData
|
||||
self.iren = iren
|
||||
self.timerId = None
|
||||
self.step = 0
|
||||
|
||||
def execute(self, obj, event):
|
||||
|
||||
print(self.timer_count)
|
||||
|
||||
dims = self.imageData.GetDimensions()
|
||||
|
||||
t=self.step/10.0
|
||||
|
||||
z0 = 2
|
||||
y0 = 4
|
||||
x0 = 4
|
||||
z1 = 14
|
||||
y1 = 12
|
||||
x1 = 12
|
||||
for z in range(dims[2]):
|
||||
for y in range(dims[1]):
|
||||
for x in range(dims[0]):
|
||||
self.imageData.SetScalarComponentFromDouble(x, y, z, 0,
|
||||
math.sin(t)*math.exp(-0.25*((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0)))
|
||||
- math.cos(t)*math.exp(-0.25*((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1))))
|
||||
|
||||
self.imageData.Modified()
|
||||
|
||||
iren = obj
|
||||
iren.GetRenderWindow().Render()
|
||||
self.timer_count += 1
|
||||
self.step += 1
|
||||
|
||||
if self.step >= self.steps :
|
||||
iren.DestroyTimer(self.timerId)
|
||||
|
||||
def WriteImage(fileName, renWin):
|
||||
'''
|
||||
'''
|
||||
|
||||
import os
|
||||
|
||||
if fileName:
|
||||
# Select the writer to use.
|
||||
path, ext = os.path.splitext(fileName)
|
||||
ext = ext.lower()
|
||||
if not ext:
|
||||
ext = '.png'
|
||||
fileName = fileName + ext
|
||||
elif ext == '.jpg':
|
||||
writer = vtkJPEGWriter()
|
||||
else:
|
||||
writer = vtkPNGWriter()
|
||||
|
||||
windowto_image_filter = vtkWindowToImageFilter()
|
||||
windowto_image_filter.SetInput(renWin)
|
||||
windowto_image_filter.SetScale(1) # image quality
|
||||
windowto_image_filter.SetInputBufferTypeToRGBA()
|
||||
|
||||
writer.SetFileName(fileName)
|
||||
writer.SetInputConnection(windowto_image_filter.GetOutputPort())
|
||||
writer.Write()
|
||||
else:
|
||||
raise RuntimeError('Need a filename.')
|
||||
|
||||
|
||||
def main():
|
||||
|
||||
colors = vtkNamedColors()
|
||||
|
||||
file_name = get_program_parameters()
|
||||
|
||||
colors.SetColor('InstantonColor', [240, 184, 160, 255])
|
||||
colors.SetColor('BackfaceColor', [255, 229, 200, 255])
|
||||
colors.SetColor('BkgColor', [51, 77, 102, 255])
|
||||
|
||||
# Create the renderer, the render window, and the interactor. The renderer
|
||||
# draws into the render window, the interactor enables mouse- and
|
||||
# keyboard-based interaction with the data within the render window.
|
||||
#
|
||||
a_renderer = vtkRenderer()
|
||||
ren_win = vtkRenderWindow()
|
||||
ren_win.AddRenderer(a_renderer)
|
||||
|
||||
iren = vtkRenderWindowInteractor()
|
||||
iren.SetRenderWindow(ren_win)
|
||||
|
||||
# The following reader is used to read a series of 2D slices (images)
|
||||
# that compose the volume. The slice dimensions are set, and the
|
||||
# pixel spacing. The data Endianness must also be specified. The reader
|
||||
# uses the FilePrefix in combination with the slice number to construct
|
||||
# filenames using the format FilePrefix.%d. (In this case the FilePrefix
|
||||
# is the root name of the file: quarter.)
|
||||
imageData = vtkImageData()
|
||||
imageData.SetDimensions(16, 16, 16)
|
||||
imageData.AllocateScalars(VTK_DOUBLE, 1)
|
||||
|
||||
dims = imageData.GetDimensions()
|
||||
|
||||
# Fill every entry of the image data with '2.0'
|
||||
# Set the input data
|
||||
for z in range(dims[2]):
|
||||
z0 = dims[2]/2
|
||||
for y in range(dims[1]):
|
||||
y0 = dims[1]/2
|
||||
for x in range(dims[0]):
|
||||
x0 = dims[0]/2
|
||||
imageData.SetScalarComponentFromDouble(x, y, z, 0, math.exp(-0.25*((x-x0)*(x-x0)+(y-y0)*(y-y0)+z*z)) - math.exp(-0.25*((x-x0)*(x-x0)+y*y+(z-z0)*(z-z0))))
|
||||
|
||||
instanton_extractor = vtkMarchingCubes()
|
||||
instanton_extractor.SetInputData(imageData)
|
||||
instanton_extractor.SetValue(0, 0.1)
|
||||
|
||||
instanton_stripper = vtkStripper()
|
||||
instanton_stripper.SetInputConnection(instanton_extractor.GetOutputPort())
|
||||
|
||||
instanton_mapper = vtkPolyDataMapper()
|
||||
instanton_mapper.SetInputConnection(instanton_stripper.GetOutputPort())
|
||||
instanton_mapper.ScalarVisibilityOff()
|
||||
|
||||
instanton = vtkActor()
|
||||
instanton.SetMapper(instanton_mapper)
|
||||
instanton.GetProperty().SetDiffuseColor(colors.GetColor3d('InstantonColor'))
|
||||
instanton.GetProperty().SetSpecular(0.3)
|
||||
instanton.GetProperty().SetSpecularPower(20)
|
||||
instanton.GetProperty().SetOpacity(0.5)
|
||||
|
||||
# The triangle stripper is used to create triangle strips from the
|
||||
# isosurface these render much faster on may systems.
|
||||
antiinstanton_extractor = vtkMarchingCubes()
|
||||
antiinstanton_extractor.SetInputData(imageData)
|
||||
antiinstanton_extractor.SetValue(0, -0.1)
|
||||
|
||||
antiinstanton_stripper = vtkStripper()
|
||||
antiinstanton_stripper.SetInputConnection(antiinstanton_extractor.GetOutputPort())
|
||||
|
||||
antiinstanton_mapper = vtkPolyDataMapper()
|
||||
antiinstanton_mapper.SetInputConnection(antiinstanton_stripper.GetOutputPort())
|
||||
antiinstanton_mapper.ScalarVisibilityOff()
|
||||
|
||||
antiinstanton = vtkActor()
|
||||
antiinstanton.SetMapper(antiinstanton_mapper)
|
||||
antiinstanton.GetProperty().SetDiffuseColor(colors.GetColor3d('Ivory'))
|
||||
|
||||
# An outline provides box around the data.
|
||||
outline_data = vtkOutlineFilter()
|
||||
outline_data.SetInputData(imageData)
|
||||
|
||||
map_outline = vtkPolyDataMapper()
|
||||
map_outline.SetInputConnection(outline_data.GetOutputPort())
|
||||
|
||||
outline = vtkActor()
|
||||
outline.SetMapper(map_outline)
|
||||
outline.GetProperty().SetColor(colors.GetColor3d('Black'))
|
||||
|
||||
# It is convenient to create an initial view of the data. The FocalPoint
|
||||
# and Position form a vector direction. Later on (ResetCamera() method)
|
||||
# this vector is used to position the camera to look at the data in
|
||||
# this direction.
|
||||
a_camera = vtkCamera()
|
||||
a_camera.SetViewUp(0, 0, -1)
|
||||
a_camera.SetPosition(0, -100, 0)
|
||||
a_camera.SetFocalPoint(0, 0, 0)
|
||||
a_camera.ComputeViewPlaneNormal()
|
||||
a_camera.Azimuth(30.0)
|
||||
a_camera.Elevation(30.0)
|
||||
|
||||
# Actors are added to the renderer. An initial camera view is created.
|
||||
# The Dolly() method moves the camera towards the FocalPoint,
|
||||
# thereby enlarging the image.
|
||||
a_renderer.AddActor(outline)
|
||||
a_renderer.AddActor(instanton)
|
||||
a_renderer.AddActor(antiinstanton)
|
||||
a_renderer.SetActiveCamera(a_camera)
|
||||
a_renderer.ResetCamera()
|
||||
a_camera.Dolly(1.0)
|
||||
|
||||
# Set a background color for the renderer and set the size of the
|
||||
# render window (expressed in pixels).
|
||||
a_renderer.SetBackground(colors.GetColor3d('BkgColor'))
|
||||
ren_win.SetSize(1024, 1024)
|
||||
ren_win.SetWindowName('ExpoDemo')
|
||||
|
||||
# Note that when camera movement occurs (as it does in the Dolly()
|
||||
# method), the clipping planes often need adjusting. Clipping planes
|
||||
# consist of two planes: near and far along the view direction. The
|
||||
# near plane clips out objects in front of the plane the far plane
|
||||
# clips out objects behind the plane. This way only what is drawn
|
||||
# between the planes is actually rendered.
|
||||
a_renderer.ResetCameraClippingRange()
|
||||
|
||||
# write image
|
||||
# WriteImage('exp.jpg',ren_win)
|
||||
|
||||
# Sign up to receive TimerEvent
|
||||
cb = vtkTimerCallback(200, imageData, iren)
|
||||
iren.AddObserver('TimerEvent', cb.execute)
|
||||
cb.timerId = iren.CreateRepeatingTimer(50)
|
||||
|
||||
# start the interaction and timer
|
||||
ren_win.Render()
|
||||
|
||||
# Initialize the event loop and then start it.
|
||||
iren.Initialize()
|
||||
iren.Start()
|
||||
|
||||
|
||||
def get_program_parameters():
|
||||
import argparse
|
||||
description = 'Simple lattice volumetric demo'
|
||||
epilogue = '''
|
||||
Derived from VTK/Examples/Cxx/Medical2.cxx
|
||||
'''
|
||||
parser = argparse.ArgumentParser(description=description, epilog=epilogue,
|
||||
formatter_class=argparse.RawDescriptionHelpFormatter)
|
||||
parser.add_argument('filename', help='FieldDensity.py')
|
||||
args = parser.parse_args()
|
||||
return args.filename
|
||||
|
||||
|
||||
def vtk_version_ok(major, minor, build):
|
||||
"""
|
||||
Check the VTK version.
|
||||
|
||||
:param major: Major version.
|
||||
:param minor: Minor version.
|
||||
:param build: Build version.
|
||||
:return: True if the requested VTK version is greater or equal to the actual VTK version.
|
||||
"""
|
||||
needed_version = 10000000000 * int(major) + 100000000 * int(minor) + int(build)
|
||||
try:
|
||||
vtk_version_number = VTK_VERSION_NUMBER
|
||||
except AttributeError: # as error:
|
||||
ver = vtkVersion()
|
||||
vtk_version_number = 10000000000 * ver.GetVTKMajorVersion() + 100000000 * ver.GetVTKMinorVersion() \
|
||||
+ ver.GetVTKBuildVersion()
|
||||
if vtk_version_number >= needed_version:
|
||||
return True
|
||||
else:
|
||||
return False
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
490
visualisation/FieldDensityAnimate.cxx
Normal file
490
visualisation/FieldDensityAnimate.cxx
Normal file
@ -0,0 +1,490 @@
|
||||
// Derived from VTK/Examples/Cxx/Medical2.cxx
|
||||
// The example reads a volume dataset, extracts two isosurfaces that
|
||||
// represent the skin and bone, and then displays them.
|
||||
//
|
||||
// Modified heavily by Peter Boyle to display lattice field theory data as movies and compare multiple files
|
||||
|
||||
#include <vtkActor.h>
|
||||
#include <vtkCamera.h>
|
||||
#include <vtkMetaImageReader.h>
|
||||
#include <vtkNamedColors.h>
|
||||
#include <vtkNew.h>
|
||||
#include <vtkOutlineFilter.h>
|
||||
#include <vtkPolyDataMapper.h>
|
||||
#include <vtkProperty.h>
|
||||
#include <vtkRenderWindow.h>
|
||||
#include <vtkRenderWindowInteractor.h>
|
||||
#include <vtkRenderer.h>
|
||||
#include <vtkStripper.h>
|
||||
#include <vtkImageData.h>
|
||||
#include <vtkVersion.h>
|
||||
#include <vtkCallbackCommand.h>
|
||||
#include <vtkTextActor.h>
|
||||
#include <vtkTextProperty.h>
|
||||
|
||||
#define MPEG
|
||||
#ifdef MPEG
|
||||
#include <vtkFFMPEGWriter.h>
|
||||
#endif
|
||||
|
||||
#include <vtkProperty2D.h>
|
||||
#include <vtkSliderWidget.h>
|
||||
#include <vtkSliderRepresentation2D.h>
|
||||
#include <vtkWindowToImageFilter.h>
|
||||
|
||||
#include <array>
|
||||
#include <string>
|
||||
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
#define USE_FLYING_EDGES
|
||||
#ifdef USE_FLYING_EDGES
|
||||
#include <vtkFlyingEdges3D.h>
|
||||
typedef vtkFlyingEdges3D isosurface;
|
||||
#else
|
||||
#include <vtkMarchingCubes.h>
|
||||
typedef vtkMarchingCubes isosurface;
|
||||
#endif
|
||||
|
||||
int mpeg = 0 ;
|
||||
int xlate = 0 ;
|
||||
|
||||
template <class T> void readFile(T& out, std::string const fname){
|
||||
#ifdef HAVE_LIME
|
||||
Grid::emptyUserRecord record;
|
||||
Grid::ScidacReader RD;
|
||||
RD.open(fname);
|
||||
RD.readScidacFieldRecord(out,record);
|
||||
RD.close();
|
||||
#endif
|
||||
}
|
||||
using namespace Grid;
|
||||
|
||||
class FrameUpdater : public vtkCallbackCommand
|
||||
{
|
||||
public:
|
||||
|
||||
FrameUpdater() {
|
||||
TimerCount = 0;
|
||||
xoff = 0;
|
||||
t = 0;
|
||||
imageData = nullptr;
|
||||
grid_data = nullptr;
|
||||
timerId = 0;
|
||||
maxCount = -1;
|
||||
}
|
||||
|
||||
static FrameUpdater* New()
|
||||
{
|
||||
FrameUpdater* cb = new FrameUpdater;
|
||||
cb->TimerCount = 0;
|
||||
return cb;
|
||||
}
|
||||
|
||||
virtual void Execute(vtkObject* caller, unsigned long eventId,void* vtkNotUsed(callData))
|
||||
{
|
||||
const int max=256;
|
||||
char text_string[max];
|
||||
|
||||
if (this->TimerCount < this->maxCount) {
|
||||
|
||||
if (vtkCommand::TimerEvent == eventId)
|
||||
{
|
||||
++this->TimerCount;
|
||||
|
||||
// Make a new frame
|
||||
auto latt_size = grid_data->Grid()->GlobalDimensions();
|
||||
for(int xx=0;xx<latt_size[0];xx++){
|
||||
for(int yy=0;yy<latt_size[1];yy++){
|
||||
for(int zz=0;zz<latt_size[2];zz++){
|
||||
int x = (xx+xoff)%latt_size[0];
|
||||
Coordinate site({x,yy,zz,t});
|
||||
RealD value = real(peekSite(*grid_data,site));
|
||||
imageData->SetScalarComponentFromDouble(xx,yy,zz,0,value);
|
||||
}}}
|
||||
|
||||
if ( xlate ) {
|
||||
xoff = (xoff + 1)%latt_size[0];
|
||||
if ( xoff== 0 ) t = (t+1)%latt_size[3];
|
||||
} else {
|
||||
t = (t+1)%latt_size[3];
|
||||
if ( t== 0 ) xoff = (xoff + 1)%latt_size[0];
|
||||
}
|
||||
|
||||
snprintf(text_string,max,"T=%d",t);
|
||||
text->SetInput(text_string);
|
||||
|
||||
std::cout << this->TimerCount<<"/"<<maxCount<< " xoff "<<xoff<<" t "<<t <<std::endl;
|
||||
imageData->Modified();
|
||||
|
||||
vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
|
||||
iren->GetRenderWindow()->Render();
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if (this->TimerCount >= this->maxCount) {
|
||||
vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
|
||||
if (this->timerId > -1)
|
||||
{
|
||||
iren->DestroyTimer(this->timerId);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
int TimerCount;
|
||||
int xoff;
|
||||
int t;
|
||||
public:
|
||||
Grid::LatticeComplexD * grid_data;
|
||||
vtkImageData* imageData = nullptr;
|
||||
vtkTextActor* text = nullptr;
|
||||
vtkFFMPEGWriter *writer = nullptr;
|
||||
int timerId ;
|
||||
int maxCount ;
|
||||
double rms;
|
||||
isosurface * posExtractor;
|
||||
isosurface * negExtractor;
|
||||
};
|
||||
class SliderCallback : public vtkCommand
|
||||
{
|
||||
public:
|
||||
static SliderCallback* New()
|
||||
{
|
||||
return new SliderCallback;
|
||||
}
|
||||
virtual void Execute(vtkObject* caller, unsigned long eventId, void* callData)
|
||||
{
|
||||
vtkSliderWidget *sliderWidget = vtkSliderWidget::SafeDownCast(caller);
|
||||
if (sliderWidget)
|
||||
{
|
||||
contour = ((vtkSliderRepresentation *)sliderWidget->GetRepresentation())->GetValue();
|
||||
}
|
||||
for(int i=0;i<fu_list.size();i++){
|
||||
fu_list[i]->posExtractor->SetValue(0, SliderCallback::contour*fu_list[i]->rms);
|
||||
fu_list[i]->negExtractor->SetValue(0, -SliderCallback::contour*fu_list[i]->rms);
|
||||
fu_list[i]->posExtractor->Modified();
|
||||
fu_list[i]->negExtractor->Modified();
|
||||
}
|
||||
}
|
||||
public:
|
||||
static double contour;
|
||||
std::vector<FrameUpdater *> fu_list;
|
||||
};
|
||||
|
||||
|
||||
double SliderCallback::contour;
|
||||
|
||||
int main(int argc, char* argv[])
|
||||
{
|
||||
using namespace Grid;
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
GridLogLayout();
|
||||
|
||||
auto latt_size = GridDefaultLatt();
|
||||
auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
|
||||
auto mpi_layout = GridDefaultMpi();
|
||||
GridCartesian Grid(latt_size, simd_layout, mpi_layout);
|
||||
|
||||
|
||||
std::cout << argc << " command Line arguments "<<std::endl;
|
||||
for(int c=0;c<argc;c++) {
|
||||
std::cout << " - "<<argv[c]<<std::endl;
|
||||
}
|
||||
|
||||
std::vector<std::string> file_list;
|
||||
double default_contour = 1.0;
|
||||
std::string arg;
|
||||
#ifdef MPEG
|
||||
if( GridCmdOptionExists(argv,argv+argc,"--mpeg") ){
|
||||
mpeg = 1;
|
||||
}
|
||||
#endif
|
||||
|
||||
if( GridCmdOptionExists(argv,argv+argc,"--xlate") ){
|
||||
xlate = 1;
|
||||
}
|
||||
|
||||
if( GridCmdOptionExists(argv,argv+argc,"--isosurface") ){
|
||||
arg=GridCmdOptionPayload(argv,argv+argc,"--isosurface");
|
||||
GridCmdOptionFloat(arg,default_contour);
|
||||
}
|
||||
if( GridCmdOptionExists(argv,argv+argc,"--file1") ){
|
||||
arg = GridCmdOptionPayload(argv,argv+argc,"--file1");
|
||||
file_list.push_back(arg);
|
||||
}
|
||||
if( GridCmdOptionExists(argv,argv+argc,"--file2") ){
|
||||
arg = GridCmdOptionPayload(argv,argv+argc,"--file2");
|
||||
file_list.push_back(arg);
|
||||
}
|
||||
if( GridCmdOptionExists(argv,argv+argc,"--file3") ){
|
||||
arg = GridCmdOptionPayload(argv,argv+argc,"--file3");
|
||||
file_list.push_back(arg);
|
||||
}
|
||||
if( GridCmdOptionExists(argv,argv+argc,"--file4") ){
|
||||
arg = GridCmdOptionPayload(argv,argv+argc,"--file4");
|
||||
file_list.push_back(arg);
|
||||
}
|
||||
for(int c=0;c<file_list.size();c++) {
|
||||
std::cout << " file: "<<file_list[c]<<std::endl;
|
||||
}
|
||||
|
||||
// Common things:
|
||||
vtkNew<vtkNamedColors> colors;
|
||||
std::array<unsigned char, 4> posColor{{240, 184, 160, 255}}; colors->SetColor("posColor", posColor.data());
|
||||
std::array<unsigned char, 4> bkg{{51, 77, 102, 255}}; colors->SetColor("BkgColor", bkg.data());
|
||||
|
||||
// Create the renderer, the render window, and the interactor. The renderer
|
||||
// draws into the render window, the interactor enables mouse- and
|
||||
// keyboard-based interaction with the data within the render window.
|
||||
//
|
||||
vtkNew<vtkRenderWindow> renWin;
|
||||
vtkNew<vtkRenderWindowInteractor> iren;
|
||||
iren->SetRenderWindow(renWin);
|
||||
|
||||
|
||||
std::vector<LatticeComplexD> data(file_list.size(),&Grid);
|
||||
FieldMetaData header;
|
||||
|
||||
|
||||
int frameCount;
|
||||
if ( mpeg ) frameCount = latt_size[3];
|
||||
else frameCount = latt_size[0] * latt_size[3];
|
||||
|
||||
std::vector<FrameUpdater *> fu_list;
|
||||
for (int f=0;f<file_list.size();f++){
|
||||
|
||||
// It is convenient to create an initial view of the data. The FocalPoint
|
||||
// and Position form a vector direction. Later on (ResetCamera() method)
|
||||
// this vector is used to position the camera to look at the data in
|
||||
// this direction.
|
||||
vtkNew<vtkCamera> aCamera;
|
||||
aCamera->SetViewUp(0, 0, -1);
|
||||
aCamera->SetPosition(0, -1000, 0);
|
||||
aCamera->SetFocalPoint(0, 0, 0);
|
||||
aCamera->ComputeViewPlaneNormal();
|
||||
aCamera->Azimuth(30.0);
|
||||
aCamera->Elevation(30.0);
|
||||
|
||||
|
||||
vtkNew<vtkRenderer> aRenderer;
|
||||
renWin->AddRenderer(aRenderer);
|
||||
|
||||
double vol = data[f].Grid()->gSites();
|
||||
|
||||
std::cout << "Reading "<<file_list[f]<<std::endl;
|
||||
readFile(data[f],file_list[f]);
|
||||
|
||||
auto nrm = norm2(data[f]);
|
||||
auto nrmbar = nrm/vol;
|
||||
auto rms = sqrt(nrmbar);
|
||||
|
||||
double contour = default_contour * rms; // default to 1 x RMS
|
||||
|
||||
// The following reader is used to read a series of 2D slices (images)
|
||||
// that compose the volume. The slice dimensions are set, and the
|
||||
// pixel spacing. The data Endianness must also be specified. The reader
|
||||
// uses the FilePrefix in combination with the slice number to construct
|
||||
// filenames using the format FilePrefix.%d. (In this case the FilePrefix
|
||||
// is the root name of the file: quarter.)
|
||||
vtkNew<vtkImageData> imageData;
|
||||
imageData->SetDimensions(latt_size[0],latt_size[1],latt_size[2]);
|
||||
imageData->AllocateScalars(VTK_DOUBLE, 1);
|
||||
for(int xx=0;xx<latt_size[0];xx++){
|
||||
for(int yy=0;yy<latt_size[1];yy++){
|
||||
for(int zz=0;zz<latt_size[2];zz++){
|
||||
Coordinate site({xx,yy,zz,0});
|
||||
RealD value = real(peekSite(data[f],site));
|
||||
imageData->SetScalarComponentFromDouble(xx,yy,zz,0,value);
|
||||
}}}
|
||||
|
||||
vtkNew<isosurface> posExtractor;
|
||||
posExtractor->SetInputData(imageData);
|
||||
posExtractor->SetValue(0, contour);
|
||||
|
||||
vtkNew<vtkStripper> posStripper;
|
||||
posStripper->SetInputConnection(posExtractor->GetOutputPort());
|
||||
|
||||
vtkNew<vtkPolyDataMapper> posMapper;
|
||||
posMapper->SetInputConnection(posStripper->GetOutputPort());
|
||||
posMapper->ScalarVisibilityOff();
|
||||
|
||||
vtkNew<vtkActor> pos;
|
||||
pos->SetMapper(posMapper);
|
||||
pos->GetProperty()->SetDiffuseColor(colors->GetColor3d("posColor").GetData());
|
||||
pos->GetProperty()->SetSpecular(0.3);
|
||||
pos->GetProperty()->SetSpecularPower(20);
|
||||
pos->GetProperty()->SetOpacity(0.5);
|
||||
|
||||
// An isosurface, or contour value is set
|
||||
// The triangle stripper is used to create triangle strips from the
|
||||
// isosurface; these render much faster on may systems.
|
||||
vtkNew<isosurface> negExtractor;
|
||||
negExtractor->SetInputData(imageData);
|
||||
negExtractor->SetValue(0, -contour);
|
||||
|
||||
vtkNew<vtkStripper> negStripper;
|
||||
negStripper->SetInputConnection(negExtractor->GetOutputPort());
|
||||
|
||||
vtkNew<vtkPolyDataMapper> negMapper;
|
||||
negMapper->SetInputConnection(negStripper->GetOutputPort());
|
||||
negMapper->ScalarVisibilityOff();
|
||||
|
||||
vtkNew<vtkActor> neg;
|
||||
neg->SetMapper(negMapper);
|
||||
neg->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());
|
||||
|
||||
// An outline provides context around the data.
|
||||
vtkNew<vtkOutlineFilter> outlineData;
|
||||
outlineData->SetInputData(imageData);
|
||||
|
||||
vtkNew<vtkPolyDataMapper> mapOutline;
|
||||
mapOutline->SetInputConnection(outlineData->GetOutputPort());
|
||||
|
||||
vtkNew<vtkActor> outline;
|
||||
outline->SetMapper(mapOutline);
|
||||
outline->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
|
||||
|
||||
vtkNew<vtkTextActor> Text;
|
||||
Text->SetInput(file_list[f].c_str());
|
||||
Text->SetPosition2(0,0);
|
||||
Text->GetTextProperty()->SetFontSize(48);
|
||||
Text->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
|
||||
|
||||
vtkNew<vtkTextActor> TextT;
|
||||
TextT->SetInput("T=0");
|
||||
TextT->SetPosition(0,.9*1025);
|
||||
TextT->GetTextProperty()->SetFontSize(48);
|
||||
TextT->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
|
||||
|
||||
|
||||
// Actors are added to the renderer. An initial camera view is created.
|
||||
// The Dolly() method moves the camera towards the FocalPoint,
|
||||
// thereby enlarging the image.
|
||||
aRenderer->AddActor(Text);
|
||||
aRenderer->AddActor(TextT);
|
||||
aRenderer->AddActor(outline);
|
||||
aRenderer->AddActor(pos);
|
||||
aRenderer->AddActor(neg);
|
||||
|
||||
// Sign up to receive TimerEvent
|
||||
vtkNew<FrameUpdater> fu;
|
||||
fu->imageData = imageData;
|
||||
fu->grid_data = &data[f];
|
||||
fu->text = TextT;
|
||||
fu->maxCount = frameCount;
|
||||
fu->posExtractor = posExtractor;
|
||||
fu->negExtractor = negExtractor;
|
||||
fu->rms = rms;
|
||||
|
||||
iren->AddObserver(vtkCommand::TimerEvent, fu);
|
||||
|
||||
aRenderer->SetActiveCamera(aCamera);
|
||||
aRenderer->ResetCamera();
|
||||
aRenderer->SetBackground(colors->GetColor3d("BkgColor").GetData());
|
||||
aCamera->Dolly(1.0);
|
||||
|
||||
double nf = file_list.size();
|
||||
std::cout << " Adding renderer " <<f<<" of "<<nf<<std::endl;
|
||||
aRenderer->SetViewport((1.0/nf)*f, 0.0,(1.0/nf)*(f+1) , 1.0);
|
||||
|
||||
// Note that when camera movement occurs (as it does in the Dolly()
|
||||
// method), the clipping planes often need adjusting. Clipping planes
|
||||
// consist of two planes: near and far along the view direction. The
|
||||
// near plane clips out objects in front of the plane; the far plane
|
||||
// clips out objects behind the plane. This way only what is drawn
|
||||
// between the planes is actually rendered.
|
||||
aRenderer->ResetCameraClippingRange();
|
||||
|
||||
fu_list.push_back(fu);
|
||||
}
|
||||
|
||||
|
||||
// Set a background color for the renderer and set the size of the
|
||||
// render window (expressed in pixels).
|
||||
// Initialize the event loop and then start it.
|
||||
renWin->SetSize(1024*file_list.size(), 1024);
|
||||
renWin->SetWindowName("FieldDensity");
|
||||
renWin->Render();
|
||||
|
||||
iren->Initialize();
|
||||
|
||||
if ( mpeg ) {
|
||||
#ifdef MPEG
|
||||
vtkWindowToImageFilter *imageFilter = vtkWindowToImageFilter::New();
|
||||
imageFilter->SetInput( renWin );
|
||||
imageFilter->SetInputBufferTypeToRGB();
|
||||
|
||||
vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New();
|
||||
writer->SetFileName("movie.avi");
|
||||
writer->SetRate(1);
|
||||
writer->SetInputConnection(imageFilter->GetOutputPort());
|
||||
writer->Start();
|
||||
|
||||
for(int i=0;i<fu_list[0]->maxCount;i++){
|
||||
for(int f=0;f<fu_list.size();f++){
|
||||
fu_list[f]->Execute(iren,vtkCommand::TimerEvent,nullptr);
|
||||
}
|
||||
imageFilter->Modified();
|
||||
writer->Write();
|
||||
}
|
||||
writer->End();
|
||||
writer->Delete();
|
||||
#else
|
||||
assert(-1 && "MPEG support not compiled");
|
||||
#endif
|
||||
} else {
|
||||
|
||||
// Add control of contour threshold
|
||||
// Create a slider widget
|
||||
vtkSmartPointer<vtkSliderRepresentation2D> sliderRep = vtkSmartPointer<vtkSliderRepresentation2D>::New();
|
||||
sliderRep->SetMinimumValue(0.1);
|
||||
sliderRep->SetMaximumValue(5.0);
|
||||
sliderRep->SetValue(1.0);
|
||||
sliderRep->SetTitleText("Fraction RMS");
|
||||
// Set color properties:
|
||||
|
||||
// Change the color of the knob that slides
|
||||
// sliderRep->GetSliderProperty()->SetColor(colors->GetColor3d("Green").GetData());
|
||||
sliderRep->GetTitleProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
|
||||
sliderRep->GetLabelProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
|
||||
sliderRep->GetSelectedProperty()->SetColor(colors->GetColor3d("DeepPink").GetData());
|
||||
|
||||
// Change the color of the bar
|
||||
sliderRep->GetTubeProperty()->SetColor(colors->GetColor3d("MistyRose").GetData());
|
||||
sliderRep->GetCapProperty()->SetColor(colors->GetColor3d("Yellow").GetData());
|
||||
sliderRep->SetSliderLength(0.05);
|
||||
sliderRep->SetSliderWidth(0.025);
|
||||
sliderRep->SetEndCapLength(0.02);
|
||||
|
||||
double nf = file_list.size();
|
||||
sliderRep->GetPoint1Coordinate()->SetCoordinateSystemToNormalizedDisplay();
|
||||
sliderRep->GetPoint1Coordinate()->SetValue(0.1, 0.1);
|
||||
sliderRep->GetPoint2Coordinate()->SetCoordinateSystemToNormalizedDisplay();
|
||||
sliderRep->GetPoint2Coordinate()->SetValue(0.9/nf, 0.1);
|
||||
|
||||
vtkSmartPointer<vtkSliderWidget> sliderWidget = vtkSmartPointer<vtkSliderWidget>::New();
|
||||
sliderWidget->SetInteractor(iren);
|
||||
sliderWidget->SetRepresentation(sliderRep);
|
||||
sliderWidget->SetAnimationModeToAnimate();
|
||||
sliderWidget->EnabledOn();
|
||||
|
||||
// Create the slider callback
|
||||
vtkSmartPointer<SliderCallback> slidercallback = vtkSmartPointer<SliderCallback>::New();
|
||||
slidercallback->fu_list = fu_list;
|
||||
sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback);
|
||||
|
||||
int timerId = iren->CreateRepeatingTimer(300);
|
||||
std::cout << "timerId: " << timerId << std::endl;
|
||||
|
||||
// Start the interaction and timer
|
||||
iren->Start();
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
113
visualisation/README
Normal file
113
visualisation/README
Normal file
@ -0,0 +1,113 @@
|
||||
========================================
|
||||
Visualisation of Grid / SciDAC format density fields using VTK (visualisation toolkit). Peter Boyle, 2025.
|
||||
========================================
|
||||
|
||||
Uses:
|
||||
|
||||
https://vtk.org
|
||||
|
||||
Files are, for example, those produced by
|
||||
Grid/HMC/ComputeWilsonFlow.cc
|
||||
and
|
||||
Grid/HMC/site_plaquette.cc
|
||||
|
||||
========================================
|
||||
Prerequisites:
|
||||
========================================
|
||||
|
||||
|
||||
1) Install ffmpeg-7.0.2 (developer install, includes headers and libraries).
|
||||
MacOS note: must install ffmpeg from source -- homebrew only installs binaries.
|
||||
|
||||
https://ffmpeg.org/download.html#releases
|
||||
|
||||
|
||||
Note: the latest ffmpeg (7.1.1) broke software compatibility with VTK.
|
||||
|
||||
|
||||
2) Build and install VTK-9.4.2, with FFMEG support enabled.
|
||||
|
||||
This is particularly involved on MacOS, so documented here.
|
||||
|
||||
cd VTK-9.4.2
|
||||
mkdir build
|
||||
cd build
|
||||
ccmake ..
|
||||
|
||||
|
||||
Using ccmake editor, set:
|
||||
|
||||
FFMPEG_DIR /usr/local
|
||||
|
||||
Toggle "advanced mode" (t)
|
||||
|
||||
Set:
|
||||
CMAKE_EXE_LINKER_FLAGS
|
||||
CMAKE_MODULE_LINKER_FLAGS
|
||||
CMAKE_SHARED_LINKER_FLAGS
|
||||
each to:
|
||||
|
||||
-framework Foundation -framework AudioToolbox -framework CoreAudio -liconv -lm -framework AVFoundation -framework CoreVideo -framework CoreMedia -framework CoreGraphics -framework AudioToolbox -framework OpenGL -framework OpenGL -framework VideoToolbox -framework CoreImage -framework AppKit -framework CoreFoundation -framework CoreServices -lz -lbz2 -Wl,-framework,CoreFoundation -Wl,-framework,Security -L/usr/local/lib -lavdevice -lavfilter -lavformat -lavcodec -lswresample -lswscale -lavutil
|
||||
|
||||
Set paths for each of
|
||||
FFMPEG_DIR /usr/local
|
||||
FFMPEG_avcodec_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_avcodec_LIBRARY /usr/local/lib/libavcodec.a
|
||||
FFMPEG_avdevice_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_avdevice_LIBRARY /usr/local/lib/libavdevice.a
|
||||
FFMPEG_avfilter_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_avfilter_LIBRARY /usr/local/lib/libavfilter.a
|
||||
FFMPEG_avformat_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_avformat_LIBRARY /usr/local/lib/libavformat.a
|
||||
FFMPEG_avresample_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_avresample_LIBRARY /usr/local/lib/libavresample.a
|
||||
FFMPEG_avutil_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_avutil_LIBRARY /usr/local/lib/libavutil.a
|
||||
FFMPEG_swresample_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_swresample_LIBRARY /usr/local/lib/libswresample.a
|
||||
FFMPEG_swscale_INCLUDE_DIR /usr/local/include
|
||||
FFMPEG_swscale_LIBRARY /usr/local/lib/libswscale.a
|
||||
|
||||
VTK_MODULE_ENABLE_VTK_IOFFMPEG YES
|
||||
|
||||
|
||||
VTK really should make it easier to pick up the flags required for FFMPEG linkage, especially as they are very quirky on MacOS.
|
||||
|
||||
|
||||
========================================
|
||||
Grid:
|
||||
========================================
|
||||
|
||||
3) Build and install a version of Grid
|
||||
|
||||
4) Ensure "grid-config" is in your path.
|
||||
|
||||
5) cd Grid/visualisation/
|
||||
libs=`grid-config --libs`
|
||||
ldflags=`grid-config --ldflags`
|
||||
cxxflags=`grid-config --cxxflags`
|
||||
cxx=`grid-config --cxx`
|
||||
|
||||
mkdir build
|
||||
cd build
|
||||
|
||||
LDFLAGS="$ldflags $libs " cmake .. -DCMAKE_CXX_COMPILER=$cxx -DCMAKE_CXX_FLAGS=$cxxflags
|
||||
|
||||
make
|
||||
|
||||
6) Invoke as:
|
||||
|
||||
FieldDensityAnimate --isosurface <float-from-0-to-5> --grid X.Y.Z.T --file1 SciDacDensityFile1 [--xlate] [--mpeg]
|
||||
FieldDensityAnimate --isosurface <float-from-0-to-5> --grid X.Y.Z.T --file1 SciDacDensityFile1 --file2 SciDacDensityFile2 [--xlate] [--mpeg]
|
||||
|
||||
==================================
|
||||
Extensions
|
||||
==================================
|
||||
|
||||
7) Direct calling from Grid ?:
|
||||
|
||||
Not yet implemented, but could develop sufficient interface to write a Lattice scalar field into MPEG direct from running code.
|
||||
|
||||
8) Example python code: FieldDensity.py . This is not interfaced to Grid.
|
||||
|
||||
|
BIN
visualisation/Topo-vs-flowtime.avi
Normal file
BIN
visualisation/Topo-vs-flowtime.avi
Normal file
Binary file not shown.
9
visualisation/cmake-command
Normal file
9
visualisation/cmake-command
Normal file
@ -0,0 +1,9 @@
|
||||
libs=`grid-config --libs`
|
||||
ldflags=`grid-config --ldflags`
|
||||
cxxflags=`grid-config --cxxflags`
|
||||
cxx=`grid-config --cxx`
|
||||
|
||||
mkdir build
|
||||
cd build
|
||||
|
||||
LDFLAGS="$ldflags $libs " cmake .. -DCMAKE_CXX_COMPILER=$cxx -DCMAKE_CXX_FLAGS=$cxxflags
|
Loading…
x
Reference in New Issue
Block a user