diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl
index 0af2dc63..fc914817 100644
--- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl
+++ b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unfitted_evolution.jl
@@ -42,7 +42,7 @@ function GridapTopOpt.evolve!(s::UnfittedFEEvolution,φh,vel,γ)
   NT, h, c, dΩ = params.NT, params.h, params.c, params.dΩ
   # ode_solver = s.ode_solver
   Tf = γ*NT
-  velh = FEFunction(V_φ,vel)
+  velh = FEFunction(V_φ,vel) # <- needs to be normalised by Hilbertian projection norm
 
   # This is temp as can't update γ in ode_solver yet
   ode_ls = LUSolver()
@@ -92,7 +92,7 @@ function GridapTopOpt.reinit!(s::UnfittedFEEvolution,φh,args...)
   reinit_nls = s.reinit_nls
 
   γd = 20
-  cₐ = 0.5 # <- 3 in connor's paper
+  cₐ = 3.0 # 0.5 # <- 3 in connor's paper
   ϵ = 1e-20
 
   # Tmp
diff --git a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unitted_evolution_test.jl b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unitted_evolution_test.jl
index b5709cf2..00a28f50 100644
--- a/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unitted_evolution_test.jl
+++ b/scripts/Embedded/MWEs/Unfitted_Evolve&Reinit/unitted_evolution_test.jl
@@ -24,16 +24,40 @@ reffe_scalar = ReferenceFE(lagrangian,Float64,order)
 V_φ = TestFESpace(model,reffe_scalar)
 Ut_φ = TransientTrialFESpace(V_φ,t -> (x->-1))
 
-φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ)
+α = 2.0*h
+a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ;
+vel_ext = VelocityExtension(a_hilb,V_φ,V_φ)
 
-ls_evo = UnfittedFEEvolution(Ut_φ, V_φ, dΩ, h; NT = 5, c = 0.1)
+# φh = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25,V_φ) # <- Already SDF
+φh = interpolate(x->cos(2π*x[1])*cos(2π*x[2])-0.11,V_φ)
 
-reinit!(ls_evo,φh)
-velh = interpolate(x->-1,V_φ)
+function compute_test_velh(φh)
+  geo = DiscreteGeometry(φh,model)
+  cutgeo = cut(model,geo)
+  Ωin = Triangulation(cutgeo,PHYSICAL_IN).a
+  diff_Ωin = DifferentiableTriangulation(Ωin)
+  dΩin = Measure(diff_Ωin,2*order)
+  j_in(φ) = ∫(φ)dΩin
+  dj_in = gradient(j_in,φh)
+  dj_vec_in = assemble_vector(dj_in,V_φ)
+  project!(vel_ext,dj_vec_in)
+  vel = dj_vec_in/sqrt(dj_vec_in'*vel_ext.K*dj_vec_in)
+  return FEFunction(V_φ,vel)
+end
+
+ls_evo = UnfittedFEEvolution(Ut_φ, V_φ, dΩ, h; NT = 5, c = 0.1,
+  reinit_nls = NLSolver(ftol=1e-12, iterations = 50, show_trace=true))
+reinit!(ls_evo,φh) # <- inital sdf
+
+velh = compute_test_velh(φh)
+# velh = interpolate(x->-1,V_φ)
 evolve!(ls_evo,φh,get_free_dof_values(velh),0.01)
 reinit!(ls_evo,φh)
 
+# You can compare the zero-iso curve for advected circle via
+# φh_expected_contour = interpolate(x->-sqrt((x[1]-0.5)^2+(x[2]-0.5)^2)+0.25+ls_evo.params.NT*0.01,V_φ)
+
 writevtk(
   Ω,"results/test_evolve",
-  cellfields=["φh"=>φh,"|∇φh|"=>norm ∘ ∇(φh)]
+  cellfields=["φh"=>φh,"|∇φh|"=>norm ∘ ∇(φh),"velh"=>velh]#,"φh_expected_contour"=>φh_expected_contour]
 )
\ No newline at end of file