diff --git a/examples/benchmark.jl b/examples/benchmark.jl
index a5b6de0..970ab63 100644
--- a/examples/benchmark.jl
+++ b/examples/benchmark.jl
@@ -2,6 +2,7 @@ using LULESH
 using MPI
 using Enzyme
 using LLVM
+using Printf
 
 LLVM.clopts("-memdep-block-scan-limit=70000")
 LLVM.clopts("-dse-memoryssa-walklimit=10000")
@@ -19,7 +20,7 @@ isdefined(Enzyme.API, :strictAliasing!) && Enzyme.API.strictAliasing!(true)
 isdefined(Enzyme.API, :typeWarning!) &&  Enzyme.API.typeWarning!(false)
 Enzyme.API.looseTypeAnalysis!(true)
 
-function main(nx, structured, num_iters, mpi, enzyme)
+function main(nx, structured, num_iters, mpi, enzyme, verification)
     # TODO: change default nr to 11
     nr = 1
     balance = 1
@@ -47,8 +48,15 @@ function main(nx, structured, num_iters, mpi, enzyme)
     domain = Domain(prob)
     if enzyme
         shadowDomain = Domain(prob)
-        @time Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
+        if verification
+            @time Enzyme.autodiff(doubleFrog, Duplicated(domain, shadowDomain))
+        else
+            @time Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
+        end
         domain = Domain(prob)
+    elseif verification
+        # not enzyme and verification means we need perturbed domain for finite differences
+        perturbedDomain = Domain(prob)
     end
 
     # getnodalMass = nodalMass(domain)
@@ -77,15 +85,49 @@ function main(nx, structured, num_iters, mpi, enzyme)
         end
     end
 
+    eps = 1e-6
+    if verification
+        num_iters = 1
+        if enzyme
+            println("seeding gradients")
+            for f in fieldnames(Domain)
+                prop = getproperty(shadowDomain, f)
+                if(typeof(prop) == Float64)
+                    setproperty!(shadowDomain,f,0.0)
+                    @printf("set %s to %f\n",f,sum(getproperty(shadowDomain,f)))
+                elseif(typeof(prop) == Vector{Float64})
+                    fill!(prop, 0.0)
+                    @printf("filled %s with %f\n",f,sum(getproperty(shadowDomain,f)))
+                end
+            end
+            fill!(shadowDomain.e, 1.0)
+        else
+            # finite difference initialization
+            println("seeding FD")
+            domain.e[1] += eps
+            perturbedDomain.e[1] -= eps
+        end
+    end
+
     # timestep to solution
     start = getWtime(prob.comm)
     while domain.time < domain.stoptime
         # this has been moved after computation of volume forces to hide launch latencies
         timeIncrement!(domain)
         if enzyme
-            Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
+            if verification
+                Enzyme.autodiff(doubleFrog, Duplicated(domain, shadowDomain))
+            else
+                Enzyme.autodiff(lagrangeLeapFrog, Duplicated(domain, shadowDomain))
+            end
         else
-            lagrangeLeapFrog(domain)
+            if verification
+                doubleFrog(domain)
+                timeIncrement!(perturbedDomain)
+                doubleFrog(perturbedDomain)
+            else
+                lagrangeLeapFrog(domain)
+            end
         end
 
         # checkErrors(domain, its, myRank)
@@ -101,6 +143,15 @@ function main(nx, structured, num_iters, mpi, enzyme)
     elapsed_time = getWtime(prob.comm) - start
     elapsed_timeG = comm_max(elapsed_time, prob.comm)
 
+    if verification
+        if enzyme
+            println("harvesting gradients")
+            @printf("gradient checksum: %f\n", shadowDomain.e[1])
+        else
+            @printf("FD checksum: %f\n",sum(domain.e - perturbedDomain.e) / (2*eps))
+        end
+    end
+
     if getMyRank(prob.comm) == 0
         verifyAndWriteFinalOutput(elapsed_timeG, domain, nx, getNumRanks(prob.comm))
     end
@@ -119,5 +170,6 @@ if !isinteractive()
     num_iters = args["num_iters"]
     mpi = args["mpi"]
     enzyme = args["enzyme"]
-    main(nx, structured, num_iters, mpi, enzyme)
+    verification = args["verification"]
+    main(nx, structured, num_iters, mpi, enzyme, verification)
 end
diff --git a/src/LULESH.jl b/src/LULESH.jl
index 2ffad9e..522ec07 100644
--- a/src/LULESH.jl
+++ b/src/LULESH.jl
@@ -56,7 +56,7 @@ include("comm.jl")
 include("utils.jl")
 
 export printUsage, IndexT, Domain, LuleshProblem, getMyRank, getNumRanks, getWtime,
-       lagrangeLeapFrog, comm_max, timeIncrement!, nodalMass, commRecv, MSG_COMM_SBN, verifyAndWriteFinalOutput
+       lagrangeLeapFrog, doubleFrog, comm_max, timeIncrement!, nodalMass, commRecv, MSG_COMM_SBN, verifyAndWriteFinalOutput
 export commSend, commRecv, commSBN
 
 function initMeshDecomp(comm)
diff --git a/src/domain.jl b/src/domain.jl
index 38c2e5c..d89ac98 100644
--- a/src/domain.jl
+++ b/src/domain.jl
@@ -3063,6 +3063,11 @@ function calcTimeConstraintsForElems(domain::Domain)
   calcHydroConstraintForElems(domain::Domain)
 end
 
+function doubleFrog(domain::Domain)
+  lagrangeLeapFrog(domain)
+  lagrangeLeapFrog(domain)
+  return nothing
+end
 
 function lagrangeLeapFrog(domain::Domain)
 
diff --git a/src/utils.jl b/src/utils.jl
index 258d234..8d5c7c3 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -21,6 +21,9 @@ function parse_cmd()
         "--enzyme"
             help = "Run Enzyme"
             action = :store_true
+        "--verification"
+            help = "Run gradient or finite difference verification"
+            action = :store_true
     end
 
     return parse_args(s)