diff --git a/CASE_3D_BFS U.jl b/CASE_3D_BFS U.jl
index 2d341119..e054fde8 100644
--- a/CASE_3D_BFS U.jl	
+++ b/CASE_3D_BFS U.jl	
@@ -12,7 +12,7 @@ mesh_file = "unv_sample_meshes/bfs_unv_tet_10mm.unv"
 
 @time mesh = UNV3D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [0.5, 0.0, 0.0]
 nu = 1e-3
@@ -22,8 +22,8 @@ model = Physics(
     time = Transient(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
-    energy = nothing,
-    domain = mesh_gpu
+    energy = Energy{Isothermal}(),
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
diff --git a/CASE_3D_BFS.jl b/CASE_3D_BFS.jl
index ddeadb80..4316666f 100644
--- a/CASE_3D_BFS.jl
+++ b/CASE_3D_BFS.jl
@@ -11,7 +11,7 @@ mesh_file = "unv_sample_meshes/bfs_unv_tet_10mm.unv"
 mesh_file = "unv_sample_meshes/bfs_OF_tet_meshes/5mm/polyMesh/"
 @time mesh = FOAM3D_mesh(mesh_file, scale=0.001, integer_type=Int64, float_type=Float64)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [0.5, 0.0, 0.0]
 nu = 1e-3
@@ -22,7 +22,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
     
 
diff --git a/CASE_3D_BFS_LES.jl b/CASE_3D_BFS_LES.jl
index 4ad434c4..8c996ff5 100644
--- a/CASE_3D_BFS_LES.jl
+++ b/CASE_3D_BFS_LES.jl
@@ -6,7 +6,7 @@ using CUDA
 mesh_file = "unv_sample_meshes/bfs_unv_hex_5mm.unv"
 mesh = UNV3D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 # INLET CONDITIONS 
 
@@ -24,7 +24,7 @@ model = Physics(
     time = Transient(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = LES{Smagorinsky}(),
-    energy = nothing,
+    energy = Energy{Isothermal}(),
     domain = mesh
     )
 
diff --git a/CASE_3D_box.jl b/CASE_3D_box.jl
index a8ea149d..2bfd2679 100644
--- a/CASE_3D_box.jl
+++ b/CASE_3D_box.jl
@@ -18,7 +18,7 @@ mesh_file="unv_sample_meshes/box_TET_PRISM_25_5mm.unv"
 
 @time mesh=UNV3D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [0.01,0.0,0.0]
 nu=1e-3
@@ -29,8 +29,8 @@ model = Physics(
     time = Steady(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
-    energy = nothing,
-    domain = mesh_gpu
+    energy = Energy{Isothermal}(),
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
diff --git a/CASE_3D_cylinder.jl b/CASE_3D_cylinder.jl
index 968c36a5..fc94a57e 100644
--- a/CASE_3D_cylinder.jl
+++ b/CASE_3D_cylinder.jl
@@ -12,7 +12,7 @@ mesh_file = "unv_sample_meshes/3D_cylinder_extruded_HEX_PRISM_FIXED.unv"
 mesh_file = "unv_sample_meshes/3D_cylinder_extruded_HEX_PRISM_FIXED_2mm_fine.unv"
 @time mesh = UNV3D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [0.50, 0.0, 0.0]
 velocity = [0.50, 0.0, 0.0]
@@ -24,8 +24,8 @@ model = Physics(
     time = Steady(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
-    energy = nothing,
-    domain = mesh_gpu
+    energy = Energy{Isothermal}(),
+    domain = mesh_dev
     )
 
 @assign! model momentum U ( 
diff --git a/CASE_3D_cylinder_U.jl b/CASE_3D_cylinder_U.jl
index d704ded5..1bdd5de2 100644
--- a/CASE_3D_cylinder_U.jl
+++ b/CASE_3D_cylinder_U.jl
@@ -10,7 +10,7 @@ mesh_file = "unv_sample_meshes/3D_cylinder.unv"
 mesh_file = "unv_sample_meshes/3D_cylinder_extruded_HEX_PRISM_FIXED_2mm.unv"
 @time mesh = UNV3D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [0.5, 0.0, 0.0]
 noSlip = [0.0, 0.0, 0.0]
@@ -21,8 +21,8 @@ model = Physics(
     time = Transient(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
-    energy = nothing,
-    domain = mesh_gpu
+    energy = Energy{Isothermal}(),
+    domain = mesh_dev
     )
 
 @assign! model momentum U ( 
diff --git a/CASE_OF_CRMHL_Wingbody.jl b/CASE_OF_CRMHL_Wingbody.jl
index 1a9d888c..fb03148a 100644
--- a/CASE_OF_CRMHL_Wingbody.jl
+++ b/CASE_OF_CRMHL_Wingbody.jl
@@ -6,8 +6,8 @@ using CUDA
 mesh_file = "unv_sample_meshes/OF_CRMHL_Wingbody_1v/polyMesh/"
 @time mesh = FOAM3D_mesh(mesh_file, scale=0.001, integer_type=Int64, float_type=Float64)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
-# mesh_gpu = mesh
+mesh_dev = adapt(CUDABackend(), mesh)
+# mesh_dev = mesh
 
 # # check volume calculation
 # volumes = ScalarField(mesh)
@@ -34,7 +34,7 @@ model = Physics(
     # turbulence = LES{Smagorinsky}(),
     energy = Energy{Isothermal}(),
     # domain = mesh
-    domain = mesh_gpu
+    domain = mesh_dev
     )
     
 walls = [:Fuselage, :FuselageAft, :Windshield, :WindshieldFrame, :WingLower, 
diff --git a/CASE_OF_buildings.jl b/CASE_OF_buildings.jl
index 2345ec3e..1dca7ab7 100644
--- a/CASE_OF_buildings.jl
+++ b/CASE_OF_buildings.jl
@@ -6,7 +6,7 @@ using CUDA
 mesh_file = "unv_sample_meshes/OF_buildings/polyMesh/"
 @time mesh = FOAM3D_mesh(mesh_file, integer_type=Int64, float_type=Float64)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 Umag = 1.5
 velocity = [Umag, 0.0, 0.0]
@@ -23,7 +23,7 @@ model = Physics(
     time = Steady(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmega}(),
-    energy = nothing,
+    energy = Energy{Isothermal}(),
     domain = mesh
     )
 
diff --git a/CASE_UNV_BFS_U.jl b/CASE_UNV_BFS_U.jl
index f370b32b..cb836d6d 100644
--- a/CASE_UNV_BFS_U.jl
+++ b/CASE_UNV_BFS_U.jl
@@ -8,7 +8,8 @@ using CUDA
 mesh_file = "unv_sample_meshes/backwardFacingStep_10mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
 
 velocity = [0.5, 0.0, 0.0]
 nu = 1e-3
@@ -18,14 +19,14 @@ model = Physics(
     time = Transient(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
-    energy = nothing,
-    domain = mesh_gpu
+    energy = Energy{Isothermal}(),
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
     Dirichlet(:inlet, velocity),
     Neumann(:outlet, 0.0),
-    Dirichlet(:wall, [0.0, 0.0, 0.0]),
+    Wall(:wall, [0.0, 0.0, 0.0]),
     Dirichlet(:top, [0.0, 0.0, 0.0])
 )
 
diff --git a/CASE_UNV_BFS_U_lowRe.jl b/CASE_UNV_BFS_U_lowRe.jl
index 20d03e00..df167203 100644
--- a/CASE_UNV_BFS_U_lowRe.jl
+++ b/CASE_UNV_BFS_U_lowRe.jl
@@ -8,7 +8,8 @@ mesh_file = "unv_sample_meshes/backwardFacingStep_5mm.unv"
 # mesh_file = "unv_sample_meshes/backwardFacingStep_2mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
 
 nu = 1e-3
 u_mag = 1.5
@@ -22,14 +23,14 @@ model = Physics(
     time = Transient(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmega}(),
-    energy = nothing,
-    domain = mesh_gpu
+    energy = Energy{Isothermal}(),
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
     Dirichlet(:inlet, velocity),
     Neumann(:outlet, 0.0),
-    Dirichlet(:wall, [0.0, 0.0, 0.0]),
+    Wall(:wall, [0.0, 0.0, 0.0]),
     Dirichlet(:top, [0.0, 0.0, 0.0])
 )
 
diff --git a/CASE_UNV_BFS_highRe.jl b/CASE_UNV_BFS_highRe.jl
index 85c5df87..b4641059 100644
--- a/CASE_UNV_BFS_highRe.jl
+++ b/CASE_UNV_BFS_highRe.jl
@@ -7,7 +7,7 @@ using CUDA
 mesh_file = "unv_sample_meshes/backwardFacingStep_5mm_long.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 nu = 1e-3
 u_mag = 3.5
@@ -24,7 +24,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmega}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
diff --git a/CASE_UNV_BFS_lowRe.jl b/CASE_UNV_BFS_lowRe.jl
index 232502a0..12b00b55 100644
--- a/CASE_UNV_BFS_lowRe.jl
+++ b/CASE_UNV_BFS_lowRe.jl
@@ -10,7 +10,8 @@ mesh_file = "unv_sample_meshes/backwardFacingStep_5mm.unv"
 mesh_file = "unv_sample_meshes/backwardFacingStep_2mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
 
 nu = 1e-3
 # u_mag = 3.5 # 2mm mesh
@@ -25,14 +26,14 @@ model = Physics(
     time = Steady(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmega}(),
-    energy = nothing,
-    domain = mesh_gpu
+    energy = Energy{Isothermal}(),
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
     Dirichlet(:inlet, velocity),
     Neumann(:outlet, 0.0),
-    Dirichlet(:wall, [0.0, 0.0, 0.0]),
+    Wall(:wall, [0.0, 0.0, 0.0]),
     Dirichlet(:top, [0.0, 0.0, 0.0])
 )
 
@@ -65,7 +66,7 @@ model = Physics(
 )
 
 schemes = (
-    U = set_schemes(gradient=Midpoint),
+    U = set_schemes(divergence=LUST, gradient=Midpoint),
     p = set_schemes(gradient=Midpoint),
     k = set_schemes(gradient=Midpoint),
     omega = set_schemes(gradient=Midpoint)
@@ -78,8 +79,8 @@ solvers = (
         preconditioner = Jacobi(),
         convergence = 1e-7,
         relax       = 0.7,
-        rtol = 1e-5,
-        atol = 1e-2
+        rtol = 1e-2,
+        atol = 1e-10
     ),
     p = set_solver(
         model.momentum.p;
@@ -87,26 +88,26 @@ solvers = (
         preconditioner = Jacobi(),
         convergence = 1e-7,
         relax       = 0.3,
-        rtol = 1e-5,
-        atol = 1e-3
+        rtol = 1e-3,
+        atol = 1e-10
     ),
     k = set_solver(
         model.turbulence.k;
         solver      = BicgstabSolver, # BicgstabSolver, GmresSolver
         preconditioner = Jacobi(),
         convergence = 1e-7,
-        relax       = 0.7,
-        rtol = 1e-5,
-        atol = 1e-2
+        relax       = 0.3,
+        rtol = 1e-2,
+        atol = 1e-10
     ),
     omega = set_solver(
         model.turbulence.omega;
         solver      = BicgstabSolver, # BicgstabSolver, GmresSolver
         preconditioner = Jacobi(),
         convergence = 1e-7,
-        relax       = 0.7,
-        rtol = 1e-5,
-        atol = 1e-2
+        relax       = 0.3,
+        rtol = 1e-2,
+        atol = 1e-10
     )
 )
 
@@ -130,7 +131,7 @@ initialise!(model.turbulence.nut, k_inlet/ω_inlet)
 
 Rx, Ry, Rz, Rp, model_out = run!(model, config) # 36.90k allocs
 
-Reff = stress_tensor(model.momentum.U, nu, model.turbulence.nut)
+Reff = stress_tensor(model.momentum.U, nu, model.turbulence.nut, config)
 Fp = pressure_force(:wall, model.momentum.p, 1.25)
 Fv = viscous_force(:wall, model.momentum.U, 1.25, nu, model.turbulence.nut)
 
diff --git a/CASE_UNV_BFS_transonic.jl b/CASE_UNV_BFS_transonic.jl
index 2b969d59..624ac304 100644
--- a/CASE_UNV_BFS_transonic.jl
+++ b/CASE_UNV_BFS_transonic.jl
@@ -6,7 +6,7 @@ using XCALibre
 mesh_file = "unv_sample_meshes/backwardFacingStep_5mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [600, 0.0, 0.0]
 nu = 1e-1
diff --git a/CASE_UNV_flatplate.jl b/CASE_UNV_flatplate.jl
index 6f62b52b..87632aef 100644
--- a/CASE_UNV_flatplate.jl
+++ b/CASE_UNV_flatplate.jl
@@ -8,8 +8,8 @@ using Krylov
 mesh_file = "unv_sample_meshes/flatplate_2D_laminar.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
-mesh_gpu = mesh
+# mesh_dev = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
 
 velocity = [0.2, 0.0, 0.0]
 nu = 1e-5
@@ -23,7 +23,7 @@ model = Physics(
         ),
     turbulence = RANS{Laminar}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
diff --git a/CASE_UNV_flatplate_LKE_ERCOFTAC.jl b/CASE_UNV_flatplate_LKE_ERCOFTAC.jl
index 9dcd94ec..d9051cdb 100644
--- a/CASE_UNV_flatplate_LKE_ERCOFTAC.jl
+++ b/CASE_UNV_flatplate_LKE_ERCOFTAC.jl
@@ -8,7 +8,7 @@ mesh_file = "unv_sample_meshes/flatplate_transition.unv"
 
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 # Turbulence Model
 velocity = [5.4,0,0]
@@ -26,7 +26,7 @@ model = Physics(
     time = Steady(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmegaLKE}(Tu = 0.01, walls=(:wall,)),
-    energy = nothing,
+    energy = Energy{Isothermal}(),
     domain = mesh
     )
 
diff --git a/CASE_UNV_flatplate_fixedT.jl b/CASE_UNV_flatplate_fixedT.jl
index a26ac0d9..797ea9f9 100644
--- a/CASE_UNV_flatplate_fixedT.jl
+++ b/CASE_UNV_flatplate_fixedT.jl
@@ -8,7 +8,7 @@ using Krylov
 mesh_file = "unv_sample_meshes/flatplate_2D_laminar.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [0.2, 0.0, 0.0]
 nu = 1e-4
diff --git a/CASE_UNV_flatplate_fixedT_KOmega.jl b/CASE_UNV_flatplate_fixedT_KOmega.jl
index 517fde00..b10f23ca 100644
--- a/CASE_UNV_flatplate_fixedT_KOmega.jl
+++ b/CASE_UNV_flatplate_fixedT_KOmega.jl
@@ -8,7 +8,7 @@ using Krylov
 mesh_file = "unv_sample_meshes/flatplate_2D_laminar.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [10, 0.0, 0.0]
 nu = 1e-4
diff --git a/CASE_UNV_flatplate_highRe.jl b/CASE_UNV_flatplate_highRe.jl
index 02c78f0c..59382406 100644
--- a/CASE_UNV_flatplate_highRe.jl
+++ b/CASE_UNV_flatplate_highRe.jl
@@ -6,8 +6,8 @@ using CUDA
 mesh_file = "unv_sample_meshes/flatplate_2D_highRe.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = mesh
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
+mesh_dev = adapt(CUDABackend(), mesh)
 
 velocity = [10, 0.0, 0.0]
 nu = 1e-5
@@ -20,7 +20,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmega}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
diff --git a/CASE_UNV_flatplate_lowRe.jl b/CASE_UNV_flatplate_lowRe.jl
index be33fd11..61d52d2c 100644
--- a/CASE_UNV_flatplate_lowRe.jl
+++ b/CASE_UNV_flatplate_lowRe.jl
@@ -6,8 +6,8 @@ using XCALibre
 mesh_file = "unv_sample_meshes/flatplate_2D_lowRe.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
-mesh_gpu = mesh
+# mesh_dev = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
 
 velocity = [10, 0.0, 0.0]
 nu = 1e-5
@@ -20,7 +20,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(),
     turbulence = RANS{KOmega}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
 
 println(typeof(model.fluid))
diff --git a/CASE_cylinder.jl b/CASE_cylinder.jl
index 65c0ddd2..1cada6ed 100644
--- a/CASE_cylinder.jl
+++ b/CASE_cylinder.jl
@@ -9,8 +9,8 @@ mesh_file = "unv_sample_meshes/cylinder_d10mm_5mm.unv"
 # mesh_file = "unv_sample_meshes/cylinder_d10mm_10-7.5-2mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
-mesh_gpu = mesh
+# mesh_dev = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
 
 # Inlet conditions
 
@@ -24,7 +24,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
 
 @assign! model momentum U ( 
diff --git a/CASE_cylinder_LES.jl b/CASE_cylinder_LES.jl
index 8789782e..b6d37a37 100644
--- a/CASE_cylinder_LES.jl
+++ b/CASE_cylinder_LES.jl
@@ -8,7 +8,7 @@ mesh_file = "unv_sample_meshes/cylinder_d10mm_2mm.unv"
 mesh_file = "unv_sample_meshes/cylinder_d10mm_10-7.5-2mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 # Inlet conditions
 
@@ -22,7 +22,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = LES{Smagorinsky}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
 
 @assign! model momentum U ( 
diff --git a/CASE_cylinder_M09.jl b/CASE_cylinder_M09.jl
index 7e88bae9..43fff088 100644
--- a/CASE_cylinder_M09.jl
+++ b/CASE_cylinder_M09.jl
@@ -9,7 +9,7 @@ mesh_file = "unv_sample_meshes/cylinder_d10mm_5mm.unv"
 # mesh_file = "unv_sample_meshes/cylinder_d10mm_10-7.5-2mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 # Inlet conditions
 
diff --git a/CASE_cylinder_RANS.jl b/CASE_cylinder_RANS.jl
index 2a295887..4f277e01 100644
--- a/CASE_cylinder_RANS.jl
+++ b/CASE_cylinder_RANS.jl
@@ -7,8 +7,8 @@ using CUDA
 mesh_file = "unv_sample_meshes/cylinder_d10mm_5mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = mesh
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = mesh
+mesh_dev = adapt(CUDABackend(), mesh)
 
 # INLET CONDITIONS 
 
@@ -28,7 +28,7 @@ model = Physics(
     # turbulence = RANS{KOmega}(β⁺=0.09),
     turbulence = RANS{KOmega}(),
     energy = Energy{Isothermal}(),
-    domain = mesh_gpu
+    domain = mesh_dev
     )
 
 @assign! model momentum U (
diff --git a/CASE_cylinder_URANS.jl b/CASE_cylinder_URANS.jl
index 6e39036f..b19eaed7 100644
--- a/CASE_cylinder_URANS.jl
+++ b/CASE_cylinder_URANS.jl
@@ -6,7 +6,7 @@ using XCALibre
 mesh_file = "unv_sample_meshes/cylinder_d10mm_5mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-mesh_gpu = adapt(CUDABackend(), mesh)
+mesh_dev = adapt(CUDABackend(), mesh)
 
 # INLET CONDITIONS 
 
@@ -24,7 +24,7 @@ model = Physics(
     time = Transient(),
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmega}(),
-    energy = nothing,
+    energy = Energy{Isothermal}(),
     domain = mesh
     )
 
diff --git a/CASE_cylinder_heated.jl b/CASE_cylinder_heated.jl
index 99f2d6b0..610e43e8 100644
--- a/CASE_cylinder_heated.jl
+++ b/CASE_cylinder_heated.jl
@@ -9,7 +9,7 @@ mesh_file = "unv_sample_meshes/cylinder_d10mm_5mm.unv"
 # mesh_file = "unv_sample_meshes/cylinder_d10mm_10-7.5-2mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 # Inlet conditions
 
diff --git a/CASE_cylinder_heated_unsteady.jl b/CASE_cylinder_heated_unsteady.jl
index 457becff..3951091b 100644
--- a/CASE_cylinder_heated_unsteady.jl
+++ b/CASE_cylinder_heated_unsteady.jl
@@ -9,7 +9,7 @@ mesh_file = "unv_sample_meshes/cylinder_d10mm_5mm.unv"
 # mesh_file = "unv_sample_meshes/cylinder_d10mm_10-7.5-2mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 # Inlet conditions
 
diff --git a/CASE_squarebend_compressible.jl b/CASE_squarebend_compressible.jl
index 3fb05852..3c8882b2 100644
--- a/CASE_squarebend_compressible.jl
+++ b/CASE_squarebend_compressible.jl
@@ -8,7 +8,7 @@ mesh_file = "unv_sample_meshes/OF_squareBend_laminar/constant/polyMesh/"
 mesh = FOAM3D_mesh(mesh_file, scale=1.0, integer_type=Int64, float_type=Float64)
 
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 # Inlet conditions
 
diff --git a/CASE_squarebend_weakly_compressible.jl b/CASE_squarebend_weakly_compressible.jl
index 0e5dff0b..84e940f8 100644
--- a/CASE_squarebend_weakly_compressible.jl
+++ b/CASE_squarebend_weakly_compressible.jl
@@ -8,7 +8,7 @@ mesh_file = "unv_sample_meshes/OF_squareBend_laminar/constant/polyMesh/"
 mesh = FOAM3D_mesh(mesh_file, scale=1.0, integer_type=Int64, float_type=Float64)
 
 
-# mesh_gpu = adapt(CUDABackend(), mesh)
+# mesh_dev = adapt(CUDABackend(), mesh)
 
 # Inlet conditions
 
diff --git a/testcases/compressible/2d_laminar_heated_plate/CASE_UNV_flatplate_fixedT.jl b/testcases/compressible/2d_laminar_heated_plate/CASE_UNV_flatplate_fixedT.jl
index 7e798c08..81ae258d 100644
--- a/testcases/compressible/2d_laminar_heated_plate/CASE_UNV_flatplate_fixedT.jl
+++ b/testcases/compressible/2d_laminar_heated_plate/CASE_UNV_flatplate_fixedT.jl
@@ -8,7 +8,7 @@ using Krylov
 # backwardFacingStep_2mm, backwardFacingStep_10mm
 mesh_file = "testcases/compressible/2d_laminar_heated_plate/flatplate_2D_laminar.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
-# mesh_gpu = adapt(CUDABackend(), mesh) # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh) # Uncomment this if using GPU
 
 velocity = [0.2, 0.0, 0.0]
 nu = 1e-4
@@ -27,7 +27,7 @@ model = Physics(
         ),
     turbulence = RANS{Laminar}(),
     energy = Energy{SensibleEnthalpy}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )
 
 @assign! model momentum U (
diff --git a/testcases/compressible/2d_turbulent_heated_plate/CASE_UNV_flatplate_fixedT_KOmega.jl b/testcases/compressible/2d_turbulent_heated_plate/CASE_UNV_flatplate_fixedT_KOmega.jl
index 36e77de5..93190d30 100644
--- a/testcases/compressible/2d_turbulent_heated_plate/CASE_UNV_flatplate_fixedT_KOmega.jl
+++ b/testcases/compressible/2d_turbulent_heated_plate/CASE_UNV_flatplate_fixedT_KOmega.jl
@@ -7,7 +7,7 @@ using Krylov
 
 mesh_file = "testcases/compressible/2d_turbulent_heated_plate/flatplate_2D_laminar.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
-# mesh_gpu = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
 
 # Inlet conditions
 velocity = [10, 0.0, 0.0]
@@ -30,7 +30,7 @@ model = Physics(
         ),
     turbulence = RANS{KOmega}(),
     energy = Energy{SensibleEnthalpy}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )
 
 @assign! model momentum U (
diff --git a/testcases/compressible/2d_unsteady_heated_cylinder/CASE_cylinder_heated_unsteady.jl b/testcases/compressible/2d_unsteady_heated_cylinder/CASE_cylinder_heated_unsteady.jl
index d719a835..b41e956a 100644
--- a/testcases/compressible/2d_unsteady_heated_cylinder/CASE_cylinder_heated_unsteady.jl
+++ b/testcases/compressible/2d_unsteady_heated_cylinder/CASE_cylinder_heated_unsteady.jl
@@ -5,7 +5,7 @@ using XCALibre
 
 mesh_file = "testcases/compressible/2d_unsteady_heated_cylinder/cylinder_d10mm_5mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
-# mesh_gpu = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
 
 # Inlet conditions
 velocity = [0.5, 0.0, 0.0]
@@ -29,7 +29,7 @@ model = Physics(
         ),
     turbulence = RANS{Laminar}(),
     energy = Energy{SensibleEnthalpy}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )
 
 @assign! model momentum U ( 
diff --git a/testcases/incompressible/2d_backwards_step/CASE_UNV_BFS.jl b/testcases/incompressible/2d_backwards_step/CASE_UNV_BFS.jl
index a239c262..e4bb803b 100644
--- a/testcases/incompressible/2d_backwards_step/CASE_UNV_BFS.jl
+++ b/testcases/incompressible/2d_backwards_step/CASE_UNV_BFS.jl
@@ -6,7 +6,7 @@ using XCALibre
 # backwardFacingStep_2mm, backwardFacingStep_10mm
 mesh_file = "testcases/incompressible/2d_backwards_step/backwardFacingStep_10mm.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
-# mesh_gpu = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
 
 # Inlet conditions
 velocity = [1.5, 0.0, 0.0]
@@ -18,7 +18,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
     energy = Energy{Isothermal}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )
 
 @assign! model momentum U (
diff --git a/testcases/incompressible/2d_laminar_plate/CASE_UNV_flatplate_lowRe.jl b/testcases/incompressible/2d_laminar_plate/CASE_UNV_flatplate_lowRe.jl
index 5b09c8ef..8f818d57 100644
--- a/testcases/incompressible/2d_laminar_plate/CASE_UNV_flatplate_lowRe.jl
+++ b/testcases/incompressible/2d_laminar_plate/CASE_UNV_flatplate_lowRe.jl
@@ -6,7 +6,7 @@ using XCALibre
 # backwardFacingStep_2mm, backwardFacingStep_10mm
 mesh_file = "testcases/incompressible/2d_laminar_plate/flatplate_2D_lowRe.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
-# mesh_gpu = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
 
 velocity = [10, 0.0, 0.0]
 nu = 1e-5
@@ -19,7 +19,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(),
     turbulence = RANS{KOmega}(),
     energy = Energy{Isothermal}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )
 
 @assign! model momentum U (
diff --git a/testcases/incompressible/2d_turbulent_plate/CASE_UNV_flatplate_highRe.jl b/testcases/incompressible/2d_turbulent_plate/CASE_UNV_flatplate_highRe.jl
index efdf4291..f9ad8605 100644
--- a/testcases/incompressible/2d_turbulent_plate/CASE_UNV_flatplate_highRe.jl
+++ b/testcases/incompressible/2d_turbulent_plate/CASE_UNV_flatplate_highRe.jl
@@ -6,7 +6,7 @@ using XCALibre
 # backwardFacingStep_2mm, backwardFacingStep_10mm
 mesh_file = "testcases/incompressible/2d_turbulent_plate/flatplate_2D_highRe.unv"
 mesh = UNV2D_mesh(mesh_file, scale=0.001)
-# mesh_gpu = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
 
 # Inlet conditions
 velocity = [10, 0.0, 0.0]
@@ -20,7 +20,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{KOmega}(),
     energy = Energy{Isothermal}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )
 
 @assign! model momentum U (
diff --git a/testcases/incompressible/3d_backwards_step/CASE_3D_BFS.jl b/testcases/incompressible/3d_backwards_step/CASE_3D_BFS.jl
index 0d7c096b..2039cc53 100644
--- a/testcases/incompressible/3d_backwards_step/CASE_3D_BFS.jl
+++ b/testcases/incompressible/3d_backwards_step/CASE_3D_BFS.jl
@@ -5,7 +5,7 @@ using XCALibre
 
 mesh_file = "testcases/incompressible/3d_backwards_step/bfs_unv_tet_10mm.unv"
 mesh = UNV3D_mesh(mesh_file, scale=0.001)
-# mesh_gpu = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
 
 # Inlet conditions
 velocity = [0.5, 0.0, 0.0]
@@ -17,7 +17,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu = nu),
     turbulence = RANS{Laminar}(),
     energy = Energy{Isothermal}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )
     
 @assign! model momentum U (
diff --git a/testcases/incompressible/3d_periodic_cascade/CASE_3D_cascade_periodic.jl b/testcases/incompressible/3d_periodic_cascade/CASE_3D_cascade_periodic.jl
index ae10b863..821fc1b5 100644
--- a/testcases/incompressible/3d_periodic_cascade/CASE_3D_cascade_periodic.jl
+++ b/testcases/incompressible/3d_periodic_cascade/CASE_3D_cascade_periodic.jl
@@ -8,7 +8,7 @@ mesh = UNV3D_mesh(mesh_file, scale=0.001)
 # backend = CUDABackend() # Uncomment this if using GPU
 backend = CPU() # Uncomment this if using CPU
 periodic = construct_periodic(mesh, backend, :top, :bottom)
-# mesh_gpu = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
+# mesh_dev = adapt(CUDABackend(), mesh)  # Uncomment this if using GPU
 
 velocity = [0.25, 0.0, 0.0]
 nu = 1e-3
@@ -19,7 +19,7 @@ model = Physics(
     fluid = Fluid{Incompressible}(nu=nu),
     turbulence = RANS{Laminar}(),
     energy = Energy{Isothermal}(),
-    domain = mesh # mesh_gpu  # use mesh_gpu for GPU backend
+    domain = mesh # mesh_dev  # use mesh_dev for GPU backend
     )