diff --git a/examples/solkit_mijin_thesis.py b/examples/solkit_mijin_thesis.py
index b1e35f7..cd9b9ba 100644
--- a/examples/solkit_mijin_thesis.py
+++ b/examples/solkit_mijin_thesis.py
@@ -173,7 +173,7 @@ def generatorSKThesis(**kwargs) -> rmk.RMKContext:
     n_b = BoundedExtrapolationDerivation("n_b")(ne)
     n_b.scalarHostProcess = rk.mpiContext.fluidProcs[-1]
 
-    extrap_LB_ui = DerivationClosure(BoundedExtrapolationDerivation("extrap_LB_ui",lowerBound=cs_b),ui)
+    extrap_LB_ui = DerivationClosure(BoundedExtrapolationDerivation("extrap_LB_ui",lowerBound=cs_b),cs)
     extrap_n = DerivationClosure(BoundedExtrapolationDerivation("extrap_n"),ni)
 
     boundaryFluxDeriv = extrap_LB_ui*extrap_n
@@ -245,7 +245,7 @@ def generatorSKThesis(**kwargs) -> rmk.RMKContext:
     rk.models.add(cm.standardBaseFluid(electronSpecies,ne,Ge,ue,Te,E,We,qe))
     rk.models.add(cm.standardBaseFluid(ionSpecies,ni,Gi,ui,Te,E))
     rk.models.add(cm.bohmBoundaryModel(electronSpecies,ne,Ge,ue,Te,cs,We,gamma))
-    rk.models.add(cm.bohmBoundaryModel(ionSpecies,ni,Gi,ui,Te,cs))
+    rk.models.add(cm.bohmBoundaryModel(ionSpecies,ni,Gi,cs,Te,cs))
     
     rk.models.add(cm.ampereMaxwell(E_dual,[Ge_dual,Gi_dual],[electronSpecies,ionSpecies],rk.norms))
 
@@ -377,7 +377,7 @@ def generatorSKThesis(**kwargs) -> rmk.RMKContext:
         for i,neut in enumerate(nn):
             neutDynModel.ddt[neut] += normConstDiff/sigmaCx[i] * stencils.DiffusionStencil(diffusionDeriv,diffCoeffOnDualGrid=True)(neut).rename("neutralDiff"+str(i+1))
 
-    neutDynModel.ddt[nn[0]] += stencils.BCDivStencil(ui,cs)(ni).rename("recyclingTerms")
+    neutDynModel.ddt[nn[0]] += stencils.BCDivStencil(cs,cs)(ni).rename("recyclingTerms")
 
     rk.models.add(neutDynModel)
     # ### CRM density and energy evolution
diff --git a/examples/solkit_mijin_thesis_kin.py b/examples/solkit_mijin_thesis_kin.py
index d2cbf6b..2895362 100644
--- a/examples/solkit_mijin_thesis_kin.py
+++ b/examples/solkit_mijin_thesis_kin.py
@@ -168,7 +168,7 @@ def generatorSKThesisKin(**kwargs) -> rmk.RMKContext:
     n_b = BoundedExtrapolationDerivation("n_b")(ni)
     n_b.scalarHostProcess = rk.mpiContext.fluidProcs[-1]
 
-    extrap_LB_ui = DerivationClosure(BoundedExtrapolationDerivation("extrap_LB_ui",lowerBound=cs_b),ui)
+    extrap_LB_ui = DerivationClosure(BoundedExtrapolationDerivation("extrap_LB_ui",lowerBound=cs_b),cs)
     extrap_n = DerivationClosure(BoundedExtrapolationDerivation("extrap_n"),ni)
 
     boundaryFluxDeriv = extrap_LB_ui*extrap_n
@@ -276,7 +276,7 @@ def generatorSKThesisKin(**kwargs) -> rmk.RMKContext:
     rk.models.add(heatingModel)
 
     rk.models.add(cm.standardBaseFluid(ionSpecies,ni,Gi,ui,Te,E))
-    rk.models.add(cm.bohmBoundaryModel(ionSpecies,ni,Gi,ui,Te,cs))
+    rk.models.add(cm.bohmBoundaryModel(ionSpecies,ni,Gi,cs,Te,cs))
 
     amModel = cm.ampereMaxwell(E_dual,[Gi_dual],[ionSpecies],rk.norms).rename("ampereMaxwell")
     amModel.ddt[E_dual] += cm.ampereMaxwellKineticElTerm(f,rk.norms).rename("kinElTerm")
@@ -337,7 +337,7 @@ def generatorSKThesisKin(**kwargs) -> rmk.RMKContext:
         for i,neut in enumerate(nn):
             neutDynModel.ddt[neut] += normConstDiff/sigmaCx[i] * stencils.DiffusionStencil(diffusionDeriv,diffCoeffOnDualGrid=True)(neut).rename("neutralDiff"+str(i+1))
 
-    neutDynModel.ddt[nn[0]] += stencils.BCDivStencil(ui,cs)(ni).rename("recyclingTerms")
+    neutDynModel.ddt[nn[0]] += stencils.BCDivStencil(cs,cs)(ni).rename("recyclingTerms")
 
     rk.models.add(neutDynModel)