temporarily add cross-check against Grid's implementation of log det …

…action to run on sunspot
lehner · Jun 12, 2024 · 75718cf · 75718cf
1 parent f82fa20
commit 75718cf
Show file tree

Hide file tree

Showing 4 changed files with 123 additions and 8 deletions.
diff --git a/lib/cgpt/lib/benchmarks.cc b/lib/cgpt/lib/benchmarks.cc
@@ -42,6 +42,80 @@ EXPORT(benchmarks,{
     return PyLong_FromLong(0);
   });
 
+#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
+#include <Grid/qcd/smearing/JacobianAction.h>
+using namespace Grid;
+
+EXPORT(test_grid,{
+    PyObject* _fields;
+    if (!PyArg_ParseTuple(args, "O", &_fields)) {
+      return NULL;
+    }
+
+    std::vector<cgpt_Lattice_base*> fields;
+    cgpt_basis_fill(fields,_fields);
+
+    PVector<LatticeColourMatrixD> cfields;
+    cgpt_basis_fill(cfields,fields);
+
+    GridCartesian* grid = (GridCartesian*)cfields[0].Grid();
+    LatticeGaugeFieldD Umu(grid);
+
+    LatticeGaugeFieldD UmuSm(grid);
+
+    LatticeGaugeFieldD Force(grid);
+
+    for(int mu=0;mu<4;mu++) {
+      PokeIndex<LorentzIndex>(Umu,cfields[mu],mu);
+      PokeIndex<LorentzIndex>(UmuSm,cfields[mu + 4],mu);
+      PokeIndex<LorentzIndex>(Force,cfields[mu + 8],mu);
+    }
+
+    typedef PeriodicGaugeImpl<GimplTypesD> PeriodicGimplD; 
+    std::cout << GridLogMessage << "P = " << WilsonLoops<PeriodicGimplD>::avgPlaquette(Umu) << std::endl;
+    std::cout << GridLogMessage << "Psm = " << WilsonLoops<PeriodicGimplD>::avgPlaquette(UmuSm) << std::endl;
+
+    // next perform smearing here and compare smeared fields
+    double rho = 0.124;
+    int Nstep  = 8;
+    typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
+    Smear_Stout<HMCWrapper::ImplPolicy> Stout(rho);
+    SmearedConfigurationMasked<HMCWrapper::ImplPolicy> SmearingPolicy(grid, Nstep, Stout);
+
+    //SmearingPolicy
+    SmearingPolicy.fill_smearedSet(Umu);
+    LatticeGaugeFieldD gridUmuSm = SmearingPolicy.get_smeared_conf(7);
+    std::cout << GridLogMessage << "PsmGrid = " << WilsonLoops<PeriodicGimplD>::avgPlaquette(gridUmuSm) << std::endl;
+
+    std::cout << GridLogMessage << "Test smeared gauge config:" << norm2(closure(gridUmuSm - UmuSm)) / norm2(UmuSm) << std::endl;
+
+    //std::cout << GridLogMessage << "Log Det Jac lvl 7:" << std::setprecision(15) << SmearingPolicy.logDetJacobian() << std::endl;
+
+    RealD ln_det = 0;
+    for (int ismr = 7; ismr > 0; --ismr) {
+      RealD d = SmearingPolicy.logDetJacobianLevel(SmearingPolicy.get_smeared_conf(ismr-1),ismr);
+      std::cout << GridLogMessage << "Log Det Jac from level" << ismr << " = " << std::setprecision(15) << d << std::endl;
+      ln_det+= d;
+    }
+    {
+      RealD d = SmearingPolicy.logDetJacobianLevel(*(SmearingPolicy.ThinLinks),0);
+      std::cout << GridLogMessage << "Log Det Jac from level" << 0 << " = " << std::setprecision(15) << d << std::endl;
+      ln_det += d;
+    }
+
+    std::cout << GridLogMessage << "Log Det Jac:" << std::setprecision(15) << -ln_det << std::endl;
+
+    LatticeGaugeFieldD gridForce(grid);
+    SmearingPolicy.logDetJacobianForce(gridForce);
+
+    //Dump(gridForce,"grid force");
+    //Dump(Force,"gpt force");
+
+    std::cout << GridLogMessage << "Test force:" << norm2(closure(gridForce - Force)) / norm2(Force) << std::endl;
+
+    return PyLong_FromLong(0);
+  });
+
 EXPORT(tests,{
     test_global_memory_system();
     return PyLong_FromLong(0);

diff --git a/lib/cgpt/lib/exports.h b/lib/cgpt/lib/exports.h
@@ -4,6 +4,7 @@ EXPORT_FUNCTION(rotate)
 EXPORT_FUNCTION(linear_combination)
 EXPORT_FUNCTION(bilinear_combination)
 EXPORT_FUNCTION(benchmarks)
+EXPORT_FUNCTION(test_grid)
 EXPORT_FUNCTION(tests)
 EXPORT_FUNCTION(create_block_map)
 EXPORT_FUNCTION(delete_block_map)

diff --git a/lib/gpt/qcd/fermion/register.py b/lib/gpt/qcd/fermion/register.py
@@ -11,24 +11,40 @@ def register(reg, op):
     reg.Mdiag = lambda dst, src: op.apply_unary_operator(2009, dst, src)
     reg.Dminus = lambda dst, src: op.apply_unary_operator(2010, dst, src)
     reg.DminusDag = lambda dst, src: op.apply_unary_operator(2011, dst, src)
-    reg.ImportPhysicalFermionSource = lambda dst, src: op.apply_unary_operator(2012, dst, src)
-    reg.ImportUnphysicalFermion = lambda dst, src: op.apply_unary_operator(2013, dst, src)
-    reg.ExportPhysicalFermionSolution = lambda dst, src: op.apply_unary_operator(2014, dst, src)
-    reg.ExportPhysicalFermionSource = lambda dst, src: op.apply_unary_operator(2015, dst, src)
+    reg.ImportPhysicalFermionSource = lambda dst, src: op.apply_unary_operator(
+        2012, dst, src
+    )
+    reg.ImportUnphysicalFermion = lambda dst, src: op.apply_unary_operator(
+        2013, dst, src
+    )
+    reg.ExportPhysicalFermionSolution = lambda dst, src: op.apply_unary_operator(
+        2014, dst, src
+    )
+    reg.ExportPhysicalFermionSource = lambda dst, src: op.apply_unary_operator(
+        2015, dst, src
+    )
     reg.Dhop = lambda dst, src: op.apply_unary_operator(3001, dst, src)
     reg.DhopDag = lambda dst, src: op.apply_unary_operator(4001, dst, src)
     reg.DhopEO = lambda dst, src: op.apply_unary_operator(3002, dst, src)
     reg.DhopEODag = lambda dst, src: op.apply_unary_operator(4002, dst, src)
-    reg.Mdir = lambda dst, src, dir, disp: op.apply_dirdisp_operator(5001, dst, src, dir, disp)
+    reg.Mdir = lambda dst, src, dir, disp: op.apply_dirdisp_operator(
+        5001, dst, src, dir, disp
+    )
     reg.MDeriv = lambda mat, dst, src: op.apply_deriv_operator(6001, mat, dst, src)
     reg.MDerivDag = lambda mat, dst, src: op.apply_deriv_operator(7001, mat, dst, src)
     reg.MoeDeriv = lambda mat, dst, src: op.apply_deriv_operator(6002, mat, dst, src)
     reg.MoeDerivDag = lambda mat, dst, src: op.apply_deriv_operator(7002, mat, dst, src)
     reg.MeoDeriv = lambda mat, dst, src: op.apply_deriv_operator(6003, mat, dst, src)
     reg.MeoDerivDag = lambda mat, dst, src: op.apply_deriv_operator(7003, mat, dst, src)
     reg.DhopDeriv = lambda mat, dst, src: op.apply_deriv_operator(6004, mat, dst, src)
-    reg.DhopDerivDag = lambda mat, dst, src: op.apply_deriv_operator(7004, mat, dst, src)
+    reg.DhopDerivDag = lambda mat, dst, src: op.apply_deriv_operator(
+        7004, mat, dst, src
+    )
     reg.DhopDerivEO = lambda mat, dst, src: op.apply_deriv_operator(6005, mat, dst, src)
-    reg.DhopDerivEODag = lambda mat, dst, src: op.apply_deriv_operator(7005, mat, dst, src)
+    reg.DhopDerivEODag = lambda mat, dst, src: op.apply_deriv_operator(
+        7005, mat, dst, src
+    )
     reg.DhopDerivOE = lambda mat, dst, src: op.apply_deriv_operator(6006, mat, dst, src)
-    reg.DhopDerivOEDag = lambda mat, dst, src: op.apply_deriv_operator(7006, mat, dst, src)
+    reg.DhopDerivOEDag = lambda mat, dst, src: op.apply_deriv_operator(
+        7006, mat, dst, src
+    )
diff --git a/tests/qcd/gauge.py b/tests/qcd/gauge.py
@@ -20,6 +20,30 @@
 # reference plaquette
 P = g.qcd.gauge.plaquette(U)
 
+# import cgpt, sys
+
+# smr = []
+# for mu in range(4):
+#     for cb in [g.even, g.odd]:
+#         smr.append(g.qcd.gauge.smear.local_stout(rho=0.124, dimension=mu, checkerboard=cb))
+
+# smrr = list(reversed(smr))
+# act = smrr[0].action_log_det_jacobian()
+# for s in smrr[1:]:
+#     act = act.transformed(s) + s.action_log_det_jacobian()
+
+# Usm = U
+# for s in smr:
+#     Usm = s(Usm)
+# g.message(P, g.qcd.gauge.plaquette(Usm), act(U))
+
+# cgpt.test_grid(U + Usm + [g(1j*x) for x in act.gradient(U, U)])
+
+# act.assert_gradient_error(rng, U, U, 1e-3, 1e-8)
+
+# sys.exit(0)
+
+
 # test rectangle calculation using parallel transport and copy_plan
 R_1x1, R_2x1 = g.qcd.gauge.rectangle(U, [(1, 1), (2, 1)])
 eps = abs(P - R_1x1)