diff --git a/docs/src/dirac.md b/docs/src/dirac.md
index 2f7fade..d0ecc3a 100644
--- a/docs/src/dirac.md
+++ b/docs/src/dirac.md
@@ -20,6 +20,16 @@ The workspace stores four fermion fields, namely `.sr`, `.sp`, `.sAp` and `.st`,
 for different purposes. If the representation is either `SU2fund` of `SU3fund`, an extra
 field with values in `U2alg`/`U3alg` is created to store the clover, used for the improvement.
 
+The functions using the fields allocated in [`DiracWorkspace`](@ref) are the following:
+
+- `dws.sr` : [`CG!`](@ref), [`flw_adapt`](@ref) (fermion case), [`bfl_error`](@ref)
+- `dws.st` : [`DwdagDw!`](@ref), [`bflw_step_vec!`](@ref)
+- `dws.sp` : [`CG!`](@ref), [`flw`](@ref) (fermion case), [`bflw_step!`](@ref), [`bflw_step_vec!`](@ref), [`propagator!`](@ref), [`bndpropagator!`](@ref), [`Tbndpropagator!`](@ref)
+- `dws.sAp` : [`CG!`](@ref), [`flw`](@ref) (fermion case), [`bflw_step!`](@ref), [`bflw_step_vec!`](@ref) 
+
+Note that other functions may call some of these functions, like [`flw_adapt`](@ref) depending on [`flw`](@ref), [`bflw!`](@ref) depending on [`bflw_step!`](@ref) or [`propagator!`](@ref) depending on [`CG!`](@ref). The fields used in the innermost function will also be modified by the outermost methods.
+
+
 ## Functions
 
 The functions [`Dw!`](@ref), [`g5Dw!`](@ref) and [`DwdagDw!`](@ref) are all related to the 
diff --git a/test/dirac/test_adapt_ferm.jl b/test/dirac/test_adapt_ferm.jl
new file mode 100644
index 0000000..0368057
--- /dev/null
+++ b/test/dirac/test_adapt_ferm.jl
@@ -0,0 +1,48 @@
+using LatticeGPU, Test, CUDA
+
+T = Float64
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4), BC_PERIODIC, (0,0,0,0,0,0))
+gp = GaugeParm{T}(SU3{T}, 6.1, 1.0)
+dpar = DiracParam{T}(SU3fund,1.3,0.9,(1.0,1.0,1.0,1.0),0.0,0.0)
+ymws = YMworkspace(SU3, T, lp)
+dws = DiracWorkspace(SU3fund,T,lp);
+
+randomize!(ymws.mom, lp, ymws)
+U = exp.(ymws.mom)
+
+psi = scalar_field(Spinor{4,SU3fund{T}},lp);
+pfrandomize!(psi,lp)
+
+Ucp = deepcopy(U)
+psicp = deepcopy(psi)
+# First Integrate very precisely up to t=2 (Wilson)
+println(" # Very precise integration ")
+wflw = wfl_rk3(Float64, 0.0004, 1.0E-7)
+flw(U,psi, wflw, 5000, gp,dpar, lp, ymws, dws)
+pl_exact = Eoft_plaq(U, gp, lp, ymws)
+cl_exact = Eoft_clover(U, gp, lp, ymws)
+println("  - Plaq:   ", pl_exact)
+println("  - Clover: ", cl_exact)
+Ufin = deepcopy(U)
+psifin = deepcopy(psi)
+
+
+# Now use Adaptive step size integrator:
+for tol in (1.0E-4, 1.0E-5, 1.0E-6, 1.0E-7, 1.0E-8)
+    local wflw = wfl_rk3(Float64, 0.0001, tol)
+    U .= Ucp
+    psi .= psicp
+    ns, eps = flw_adapt(U,psi, wflw, 2.0, gp,dpar,lp, ymws,dws)
+    pl = Eoft_plaq(U, gp, lp, ymws)
+    cl = Eoft_clover(U, gp, lp, ymws)
+    psierr = sum(norm2.((psi.-psifin)))./prod(lp.iL)
+
+    println(" # Adaptive integrator (tol=$tol): ", ns, " steps")
+    U .= U ./ Ufin
+    maxd = CUDA.mapreduce(dev_one, max, U, init=0.0)
+    println("  - Plaq:   ", pl," [diff: ", abs(pl-pl_exact), "; ",
+            maxd, "]")
+    println("  - Clover: ", cl, " [diff: ", abs(cl-cl_exact), "; ",
+            maxd, "]")
+    println("  - Fermion diff: ", psierr)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 31fea0b..6ffdddf 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,4 +7,4 @@ include("dirac/test_solver_rand.jl")
 include("dirac/test_flow_tl.jl")
 include("dirac/test_backflow_tl.jl")
 include("dirac/test_backflow.jl")
-
+include("dirac/test_adapt_ferm.jl")