diff --git a/src/analysis/sampling/warmup.jl b/src/analysis/sampling/warmup.jl
index 222348b091f48f6514b646182db56229e1af1213..a9768eb62ad58c1a7d93ed48e72d074e5c73670e 100644
--- a/src/analysis/sampling/warmup.jl
+++ b/src/analysis/sampling/warmup.jl
@@ -16,17 +16,22 @@ function warmup(
     modifications = [],
     workerids = [myid()],
 )
-    # create optimization problem, apply constraints, load on all workers
+    # create optimization problem, apply constraints
     save_model = :(
         begin
             optmodel = $COBREXA.make_optimization_model($model, $optimizer)
+            for mod in $modifications
+                mod($model, optmodel)
+            end
             optmodel
         end
     )
 
+    # load on all workers
     map(fetch, save_at.(workerids, :cobrexa_sampling_warmup_optmodel, Ref(save_model)))
-
-    ret = m -> value.(m[:x]) # get all the fluxes
+    
+    # generate warm up points, like FVA
+    ret = m -> value.(m[:x]) # get all the fluxes, not just max/min like FVA
     fluxes = dpmap(
         rid -> :($COBREXA._FVA_optimize_reaction(
             cobrexa_sampling_warmup_optmodel,
diff --git a/test/analysis/sampling/warmup.jl b/test/analysis/sampling/warmup.jl
new file mode 100644
index 0000000000000000000000000000000000000000..3a2e894f5d5022b81f11a7bcd6c7d3d74b6b6f19
--- /dev/null
+++ b/test/analysis/sampling/warmup.jl
@@ -0,0 +1,48 @@
+@testset "Warm up point generation" begin
+    model_path = download_data_file(
+        "http://bigg.ucsd.edu/static/models/e_coli_core.json",
+        joinpath("data", "e_coli_core.json"),
+        "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
+    )
+    model = load_model(StandardModel, model_path)
+
+    using Distributed
+
+    nps = 4 - nprocs()
+
+    # load extra processes, have at least 4 available
+    if 1 <= nps <= 3 
+        addprocs(nps)
+    end
+    @everywhere using COBREXA, Tulip
+
+    # Serial test
+    ws, lbs, ubs = warmup(
+        model,
+        Tulip.Optimizer;
+        modifications = [change_constraint("EX_glc__D_e", -4, 4)],
+        warmup_points = collect(1:n_reactions(model)),
+        workerids = [myid()],
+    )
+
+    ind = first(indexin(["EX_glc__D_e"], reactions(model)))
+    @test size(ws) == (95, 2)
+    @test size(ws[1,1]) == (95,)
+    @test lbs[ind] ≈ -4
+    @test ubs[ind] ≈ 4
+    
+    # Parallel test
+    ws, lbs, ubs = warmup(
+        model,
+        Tulip.Optimizer;
+        modifications = [change_constraint("EX_glc__D_e", -4, 4)],
+        warmup_points = collect(1:n_reactions(model)),
+        workerids = [workers()],
+    )
+
+    ind = first(indexin(["EX_glc__D_e"], reactions(model)))
+    @test size(ws) == (95, 2)
+    @test size(ws[1,1]) == (95,)
+    @test lbs[ind] ≈ -4
+    @test ubs[ind] ≈ 4
+end