mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-09 21:50:45 +01:00
Implement EOFA pseudofermion force and Shamir tests for G-parity and non G-parity cases
This commit is contained in:
David Murphy
parent
ec1e2f7a40
commit
47a12ec7b5
@ -118,7 +118,7 @@ namespace QCD{
|
|||||||
RealD scale = std::sqrt(0.5);
|
RealD scale = std::sqrt(0.5);
|
||||||
gaussian(pRNG,eta);
|
gaussian(pRNG,eta);
|
||||||
eta = eta * scale;
|
eta = eta * scale;
|
||||||
printf("Heatbath source vector: <\eta|\eta> = %1.15e\n", norm2(eta));
|
printf("Heatbath source vector: <\\eta|\\eta> = %1.15e\n", norm2(eta));
|
||||||
|
|
||||||
// \Phi = ( \alpha_{0} + \sum_{k=1}^{N_{p}} \alpha_{l} * \gamma_{l} ) * \eta
|
// \Phi = ( \alpha_{0} + \sum_{k=1}^{N_{p}} \alpha_{l} * \gamma_{l} ) * \eta
|
||||||
RealD N(PowerNegHalf.norm);
|
RealD N(PowerNegHalf.norm);
|
||||||
@ -223,6 +223,40 @@ namespace QCD{
|
|||||||
// EOFA pseudofermion force: see Eqns. (34)-(36) of arXiv:1706.05843
|
// EOFA pseudofermion force: see Eqns. (34)-(36) of arXiv:1706.05843
|
||||||
virtual void deriv(const GaugeField& U, GaugeField& dSdU)
|
virtual void deriv(const GaugeField& U, GaugeField& dSdU)
|
||||||
{
|
{
|
||||||
|
Lop.ImportGauge(U);
|
||||||
|
Rop.ImportGauge(U);
|
||||||
|
|
||||||
|
FermionField spProj_Phi (Lop.FermionGrid());
|
||||||
|
FermionField Omega_spProj_Phi(Lop.FermionGrid());
|
||||||
|
FermionField CG_src (Lop.FermionGrid());
|
||||||
|
FermionField Chi (Lop.FermionGrid());
|
||||||
|
FermionField g5_R5_Chi (Lop.FermionGrid());
|
||||||
|
|
||||||
|
GaugeField force(Lop.GaugeGrid());
|
||||||
|
|
||||||
|
// LH: dSdU = k \chi_{L}^{\dagger} \gamma_{5} R_{5} ( \partial_{x,\mu} D_{w} ) \chi_{L}
|
||||||
|
// \chi_{L} = H(mf)^{-1} \Omega_{-} P_{-} \Phi
|
||||||
|
spProj(Phi, spProj_Phi, -1, Lop.Ls);
|
||||||
|
Lop.Omega(spProj_Phi, Omega_spProj_Phi, -1, 0);
|
||||||
|
G5R5(CG_src, Omega_spProj_Phi);
|
||||||
|
spProj_Phi = zero;
|
||||||
|
Solver(Lop, CG_src, spProj_Phi);
|
||||||
|
Lop.Dtilde(spProj_Phi, Chi);
|
||||||
|
G5R5(g5_R5_Chi, Chi);
|
||||||
|
Lop.MDeriv(force, g5_R5_Chi, Chi, DaggerNo);
|
||||||
|
dSdU = Lop.k * force;
|
||||||
|
|
||||||
|
// RH: dSdU = dSdU - k \chi_{R}^{\dagger} \gamma_{5} R_{5} ( \partial_{x,\mu} D_{w} ) \chi_{}
|
||||||
|
// \chi_{R} = ( H(mb) - \Delta_{+}(mf,mb) P_{+} )^{-1} \Omega_{+} P_{+} \Phi
|
||||||
|
spProj(Phi, spProj_Phi, 1, Rop.Ls);
|
||||||
|
Rop.Omega(spProj_Phi, Omega_spProj_Phi, 1, 0);
|
||||||
|
G5R5(CG_src, Omega_spProj_Phi);
|
||||||
|
spProj_Phi = zero;
|
||||||
|
Solver(Rop, CG_src, spProj_Phi);
|
||||||
|
Rop.Dtilde(spProj_Phi, Chi);
|
||||||
|
G5R5(g5_R5_Chi, Chi);
|
||||||
|
Lop.MDeriv(force, g5_R5_Chi, Chi, DaggerNo);
|
||||||
|
dSdU = dSdU - Rop.k * force;
|
||||||
};
|
};
|
||||||
};
|
};
|
||||||
}}
|
}}
|
||||||
|
|||||||
164
tests/forces/Test_dwf_force_eofa.cc
Normal file
164
tests/forces/Test_dwf_force_eofa.cc
Normal file
@ -0,0 +1,164 @@
|
|||||||
|
/*************************************************************************************
|
||||||
|
|
||||||
|
Grid physics library, www.github.com/paboyle/Grid
|
||||||
|
|
||||||
|
Source file: ./tests/forces/Test_dwf_force_eofa.cc
|
||||||
|
|
||||||
|
Copyright (C) 2017
|
||||||
|
|
||||||
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||||
|
Author: David Murphy <dmurphy@phys.columbia.edu>
|
||||||
|
|
||||||
|
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;
|
||||||
|
using namespace Grid::QCD;
|
||||||
|
|
||||||
|
int main (int argc, char** argv)
|
||||||
|
{
|
||||||
|
Grid_init(&argc, &argv);
|
||||||
|
|
||||||
|
std::vector<int> latt_size = GridDefaultLatt();
|
||||||
|
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||||
|
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||||
|
|
||||||
|
const int Ls = 8;
|
||||||
|
|
||||||
|
GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi());
|
||||||
|
GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||||
|
GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid);
|
||||||
|
GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid);
|
||||||
|
|
||||||
|
// Want a different conf at every run
|
||||||
|
// First create an instance of an engine.
|
||||||
|
std::random_device rnd_device;
|
||||||
|
// Specify the engine and distribution.
|
||||||
|
std::mt19937 mersenne_engine(rnd_device());
|
||||||
|
std::uniform_int_distribution<int> dist(1, 100);
|
||||||
|
|
||||||
|
auto gen = std::bind(dist, mersenne_engine);
|
||||||
|
std::vector<int> seeds4(4);
|
||||||
|
generate(begin(seeds4), end(seeds4), gen);
|
||||||
|
|
||||||
|
//std::vector<int> seeds4({1,2,3,5});
|
||||||
|
std::vector<int> seeds5({5,6,7,8});
|
||||||
|
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||||
|
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||||
|
|
||||||
|
int threads = GridThread::GetThreads();
|
||||||
|
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
|
||||||
|
|
||||||
|
LatticeFermion phi (FGrid); gaussian(RNG5, phi);
|
||||||
|
LatticeFermion Mphi (FGrid);
|
||||||
|
LatticeFermion MphiPrime (FGrid);
|
||||||
|
|
||||||
|
LatticeGaugeField U(UGrid);
|
||||||
|
SU3::HotConfiguration(RNG4,U);
|
||||||
|
|
||||||
|
////////////////////////////////////
|
||||||
|
// Unmodified matrix element
|
||||||
|
////////////////////////////////////
|
||||||
|
RealD mf = 0.01;
|
||||||
|
RealD mb = 1.0;
|
||||||
|
RealD M5 = 1.8;
|
||||||
|
DomainWallEOFAFermionR Lop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, 0.0, -1, M5);
|
||||||
|
DomainWallEOFAFermionR Rop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, -1.0, 1, M5);
|
||||||
|
OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, 12);
|
||||||
|
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
|
||||||
|
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true);
|
||||||
|
|
||||||
|
Meofa.refresh(U, RNG5);
|
||||||
|
RealD S = Meofa.S(U); // pdag M p
|
||||||
|
|
||||||
|
// get the deriv of phidag M phi with respect to "U"
|
||||||
|
LatticeGaugeField UdSdU(UGrid);
|
||||||
|
Meofa.deriv(U, UdSdU);
|
||||||
|
|
||||||
|
////////////////////////////////////
|
||||||
|
// Modify the gauge field a little
|
||||||
|
////////////////////////////////////
|
||||||
|
RealD dt = 0.0001;
|
||||||
|
|
||||||
|
LatticeColourMatrix mommu(UGrid);
|
||||||
|
LatticeColourMatrix forcemu(UGrid);
|
||||||
|
LatticeGaugeField mom(UGrid);
|
||||||
|
LatticeGaugeField Uprime(UGrid);
|
||||||
|
|
||||||
|
for(int mu=0; mu<Nd; mu++){
|
||||||
|
|
||||||
|
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
|
||||||
|
|
||||||
|
PokeIndex<LorentzIndex>(mom, mommu, mu);
|
||||||
|
|
||||||
|
// fourth order exponential approx
|
||||||
|
parallel_for(auto i=mom.begin(); i<mom.end(); i++){
|
||||||
|
Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt + mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*Ddwf.ImportGauge(Uprime);
|
||||||
|
Ddwf.M (phi,MphiPrime);
|
||||||
|
|
||||||
|
ComplexD Sprime = innerProduct(MphiPrime ,MphiPrime);*/
|
||||||
|
RealD Sprime = Meofa.S(Uprime);
|
||||||
|
|
||||||
|
//////////////////////////////////////////////
|
||||||
|
// Use derivative to estimate dS
|
||||||
|
//////////////////////////////////////////////
|
||||||
|
|
||||||
|
LatticeComplex dS(UGrid);
|
||||||
|
dS = zero;
|
||||||
|
for(int mu=0; mu<Nd; mu++){
|
||||||
|
mommu = PeekIndex<LorentzIndex>(UdSdU, mu);
|
||||||
|
mommu = Ta(mommu)*2.0;
|
||||||
|
PokeIndex<LorentzIndex>(UdSdU, mommu, mu);
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int mu=0; mu<Nd; mu++){
|
||||||
|
forcemu = PeekIndex<LorentzIndex>(UdSdU, mu);
|
||||||
|
mommu = PeekIndex<LorentzIndex>(mom, mu);
|
||||||
|
|
||||||
|
// Update PF action density
|
||||||
|
dS = dS + trace(mommu*forcemu)*dt;
|
||||||
|
}
|
||||||
|
|
||||||
|
ComplexD dSpred = sum(dS);
|
||||||
|
|
||||||
|
/*std::cout << GridLogMessage << " S " << S << std::endl;
|
||||||
|
std::cout << GridLogMessage << " Sprime " << Sprime << std::endl;
|
||||||
|
std::cout << GridLogMessage << "dS " << Sprime-S << std::endl;
|
||||||
|
std::cout << GridLogMessage << "predict dS " << dSpred << std::endl;*/
|
||||||
|
printf("\nS = %1.15e\n", S);
|
||||||
|
printf("Sprime = %1.15e\n", Sprime);
|
||||||
|
printf("dS = %1.15e\n", Sprime - S);
|
||||||
|
printf("real(dS_predict) = %1.15e\n", dSpred.real());
|
||||||
|
printf("imag(dS_predict) = %1.15e\n\n", dSpred.imag());
|
||||||
|
|
||||||
|
assert( fabs(real(Sprime-S-dSpred)) < 1.0 ) ;
|
||||||
|
|
||||||
|
std::cout << GridLogMessage << "Done" << std::endl;
|
||||||
|
Grid_finalize();
|
||||||
|
}
|
||||||
169
tests/forces/Test_dwf_gpforce_eofa.cc
Normal file
169
tests/forces/Test_dwf_gpforce_eofa.cc
Normal file
@ -0,0 +1,169 @@
|
|||||||
|
/*************************************************************************************
|
||||||
|
|
||||||
|
Grid physics library, www.github.com/paboyle/Grid
|
||||||
|
|
||||||
|
Source file: ./tests/forces/Test_dwf_force_eofa.cc
|
||||||
|
|
||||||
|
Copyright (C) 2017
|
||||||
|
|
||||||
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||||
|
Author: David Murphy <dmurphy@phys.columbia.edu>
|
||||||
|
|
||||||
|
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;
|
||||||
|
using namespace Grid::QCD;
|
||||||
|
|
||||||
|
typedef GparityWilsonImplR FermionImplPolicy;
|
||||||
|
typedef GparityDomainWallEOFAFermionR FermionAction;
|
||||||
|
typedef typename FermionAction::FermionField FermionField;
|
||||||
|
|
||||||
|
int main (int argc, char** argv)
|
||||||
|
{
|
||||||
|
Grid_init(&argc, &argv);
|
||||||
|
|
||||||
|
std::vector<int> latt_size = GridDefaultLatt();
|
||||||
|
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||||
|
std::vector<int> mpi_layout = GridDefaultMpi();
|
||||||
|
|
||||||
|
const int Ls = 8;
|
||||||
|
|
||||||
|
GridCartesian *UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()), GridDefaultMpi());
|
||||||
|
GridRedBlackCartesian *UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||||
|
GridCartesian *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid);
|
||||||
|
GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid);
|
||||||
|
|
||||||
|
// Want a different conf at every run
|
||||||
|
// First create an instance of an engine.
|
||||||
|
std::random_device rnd_device;
|
||||||
|
// Specify the engine and distribution.
|
||||||
|
std::mt19937 mersenne_engine(rnd_device());
|
||||||
|
std::uniform_int_distribution<int> dist(1, 100);
|
||||||
|
|
||||||
|
auto gen = std::bind(dist, mersenne_engine);
|
||||||
|
std::vector<int> seeds4(4);
|
||||||
|
generate(begin(seeds4), end(seeds4), gen);
|
||||||
|
|
||||||
|
//std::vector<int> seeds4({1,2,3,5});
|
||||||
|
std::vector<int> seeds5({5,6,7,8});
|
||||||
|
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||||
|
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||||
|
|
||||||
|
int threads = GridThread::GetThreads();
|
||||||
|
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
|
||||||
|
|
||||||
|
LatticeFermion phi (FGrid); gaussian(RNG5, phi);
|
||||||
|
LatticeFermion Mphi (FGrid);
|
||||||
|
LatticeFermion MphiPrime (FGrid);
|
||||||
|
|
||||||
|
LatticeGaugeField U(UGrid);
|
||||||
|
SU3::HotConfiguration(RNG4,U);
|
||||||
|
|
||||||
|
////////////////////////////////////
|
||||||
|
// Unmodified matrix element
|
||||||
|
////////////////////////////////////
|
||||||
|
RealD mf = 0.01;
|
||||||
|
RealD mb = 1.0;
|
||||||
|
RealD M5 = 1.8;
|
||||||
|
FermionAction::ImplParams params;
|
||||||
|
FermionAction Lop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, 0.0, -1, M5, params);
|
||||||
|
FermionAction Rop(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, -1.0, 1, M5, params);
|
||||||
|
OneFlavourRationalParams Params(0.95, 100.0, 5000, 1.0e-12, 12);
|
||||||
|
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
|
||||||
|
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
|
||||||
|
|
||||||
|
Meofa.refresh(U, RNG5);
|
||||||
|
RealD S = Meofa.S(U); // pdag M p
|
||||||
|
|
||||||
|
// get the deriv of phidag M phi with respect to "U"
|
||||||
|
LatticeGaugeField UdSdU(UGrid);
|
||||||
|
Meofa.deriv(U, UdSdU);
|
||||||
|
|
||||||
|
////////////////////////////////////
|
||||||
|
// Modify the gauge field a little
|
||||||
|
////////////////////////////////////
|
||||||
|
RealD dt = 0.0001;
|
||||||
|
|
||||||
|
LatticeColourMatrix mommu(UGrid);
|
||||||
|
LatticeColourMatrix forcemu(UGrid);
|
||||||
|
LatticeGaugeField mom(UGrid);
|
||||||
|
LatticeGaugeField Uprime(UGrid);
|
||||||
|
|
||||||
|
for(int mu=0; mu<Nd; mu++){
|
||||||
|
|
||||||
|
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
|
||||||
|
|
||||||
|
PokeIndex<LorentzIndex>(mom, mommu, mu);
|
||||||
|
|
||||||
|
// fourth order exponential approx
|
||||||
|
parallel_for(auto i=mom.begin(); i<mom.end(); i++){
|
||||||
|
Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt + mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
|
||||||
|
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/*Ddwf.ImportGauge(Uprime);
|
||||||
|
Ddwf.M (phi,MphiPrime);
|
||||||
|
|
||||||
|
ComplexD Sprime = innerProduct(MphiPrime ,MphiPrime);*/
|
||||||
|
RealD Sprime = Meofa.S(Uprime);
|
||||||
|
|
||||||
|
//////////////////////////////////////////////
|
||||||
|
// Use derivative to estimate dS
|
||||||
|
//////////////////////////////////////////////
|
||||||
|
|
||||||
|
LatticeComplex dS(UGrid);
|
||||||
|
dS = zero;
|
||||||
|
for(int mu=0; mu<Nd; mu++){
|
||||||
|
mommu = PeekIndex<LorentzIndex>(UdSdU, mu);
|
||||||
|
mommu = Ta(mommu)*2.0;
|
||||||
|
PokeIndex<LorentzIndex>(UdSdU, mommu, mu);
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int mu=0; mu<Nd; mu++){
|
||||||
|
forcemu = PeekIndex<LorentzIndex>(UdSdU, mu);
|
||||||
|
mommu = PeekIndex<LorentzIndex>(mom, mu);
|
||||||
|
|
||||||
|
// Update PF action density
|
||||||
|
dS = dS + trace(mommu*forcemu)*dt;
|
||||||
|
}
|
||||||
|
|
||||||
|
ComplexD dSpred = sum(dS);
|
||||||
|
|
||||||
|
/*std::cout << GridLogMessage << " S " << S << std::endl;
|
||||||
|
std::cout << GridLogMessage << " Sprime " << Sprime << std::endl;
|
||||||
|
std::cout << GridLogMessage << "dS " << Sprime-S << std::endl;
|
||||||
|
std::cout << GridLogMessage << "predict dS " << dSpred << std::endl;*/
|
||||||
|
printf("\nS = %1.15e\n", S);
|
||||||
|
printf("Sprime = %1.15e\n", Sprime);
|
||||||
|
printf("dS = %1.15e\n", Sprime - S);
|
||||||
|
printf("real(dS_predict) = %1.15e\n", dSpred.real());
|
||||||
|
printf("imag(dS_predict) = %1.15e\n\n", dSpred.imag());
|
||||||
|
|
||||||
|
assert( fabs(real(Sprime-S-dSpred)) < 1.0 ) ;
|
||||||
|
|
||||||
|
std::cout << GridLogMessage << "Done" << std::endl;
|
||||||
|
Grid_finalize();
|
||||||
|
}
|
||||||
Loading…
x
Reference in New Issue
Block a user