diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 417737777cae37ddb51f595203d6aa3b261ca9b6..0236d8d7995bedb254062d6dd29f26bd064b04b1 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -15,9 +15,6 @@
 
 #include <Kokkos_Core.hpp>
 
-#include <array>
-#include <cstdio>
-
 template <typename T>
 class MeshTransformation;
 template <typename OutputType, typename... InputType>
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 74444564248a26a25c9062dcd9ce5e38c6e21ed6..4e0f8f73788bd46320e03ae6ad8a9810829efd0e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -14,8 +14,52 @@
 
 /////////// TEMPORARY
 
+#include <language/utils/PugsFunctionAdapter.hpp>
 #include <output/VTKWriter.hpp>
 
+template <typename T>
+class InterpolateItemValue;
+template <typename OutputType, typename InputType>
+class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+  static constexpr size_t Dimension = OutputType::Dimension;
+  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
+
+ public:
+  template <ItemType item_type>
+  static inline ItemValue<OutputType, item_type>
+  interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
+  {
+    const auto flatten_args = Adapter::getFlattenArgs(function_symbol_id);
+
+    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
+    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+    const IConnectivity& connectivity = *position.connectivity_ptr();
+
+    ItemValue<OutputType, item_type> value(connectivity);
+    using ItemId = ItemIdT<item_type>;
+
+    parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &flatten_args, &tokens](ItemId i) {
+      const int32_t t = tokens.acquire();
+
+      auto& execution_policy = context_list[t];
+
+      Adapter::convertArgs(execution_policy.currentContext(), flatten_args, position[i]);
+      auto result = expression.execute(execution_policy);
+      value[i]    = convert_result(std::move(result));
+
+      tokens.release(t);
+    });
+
+    return value;
+  }
+};
+
 template <size_t Dimension>
 struct GlaceScheme
 {
@@ -27,7 +71,10 @@ struct GlaceScheme
   const MeshType& m_mesh;
 
   GlaceScheme(const IMesh& mesh,
-              const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+              const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+              FunctionSymbolId rho_id,
+              FunctionSymbolId u_id,
+              FunctionSymbolId p_id)
     : m_mesh{dynamic_cast<const MeshType&>(mesh)}
   {
     MeshDataType mesh_data(m_mesh);
@@ -71,7 +118,18 @@ struct GlaceScheme
     }
 
     UnknownsType unknowns(mesh_data);
-    unknowns.initializeSod();
+
+    unknowns.rhoj() =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, mesh_data.xj());
+
+    unknowns.pj() =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
+
+    unknowns.uj() =
+      InterpolateItemValue<TinyVector<Dimension>(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id,
+                                                                                                               mesh_data
+                                                                                                                 .xj());
+    unknowns.gammaj().fill(1.4);
 
     AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
@@ -83,13 +141,27 @@ struct GlaceScheme
     int itermax   = std::numeric_limits<int>::max();
     int iteration = 0;
 
-    CellValue<double>& rhoj   = unknowns.rhoj();
-    CellValue<double>& ej     = unknowns.ej();
-    CellValue<double>& pj     = unknowns.pj();
-    CellValue<double>& gammaj = unknowns.gammaj();
-    CellValue<double>& cj     = unknowns.cj();
+    CellValue<double>& rhoj              = unknowns.rhoj();
+    CellValue<double>& ej                = unknowns.ej();
+    CellValue<double>& pj                = unknowns.pj();
+    CellValue<double>& gammaj            = unknowns.gammaj();
+    CellValue<double>& cj                = unknowns.cj();
+    CellValue<TinyVector<Dimension>>& uj = unknowns.uj();
+    CellValue<double>& Ej                = unknowns.Ej();
+    CellValue<double>& mj                = unknowns.mj();
+    CellValue<double>& inv_mj            = unknowns.invMj();
 
     BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+    block_eos.updateEandCFromRhoP();
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
 
     VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
 
@@ -170,35 +242,37 @@ SchemeModule::SchemeModule()
 
                             ));
 
-  this
-    ->_addBuiltinFunction("glace",
-                          std::make_shared<
-                            BuiltinFunctionEmbedder<void, std::shared_ptr<const IMesh>,
-                                                    std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>>>(
-                            std::function<void(std::shared_ptr<const IMesh>,
-                                               std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>)>{
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  GlaceScheme<1>{*p_mesh, bc_descriptor_list};
-                                  break;
-                                }
-                                case 2: {
-                                  GlaceScheme<2>{*p_mesh, bc_descriptor_list};
-                                  break;
-                                }
-                                case 3: {
-                                  GlaceScheme<3>{*p_mesh, bc_descriptor_list};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
+  this->_addBuiltinFunction("glace",
+                            std::make_shared<
+                              BuiltinFunctionEmbedder<void, std::shared_ptr<const IMesh>,
+                                                      std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>,
+                                                      FunctionSymbolId, FunctionSymbolId, FunctionSymbolId>>(
+                              std::function<void(std::shared_ptr<const IMesh>,
+                                                 std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>,
+                                                 FunctionSymbolId, FunctionSymbolId, FunctionSymbolId)>{
+
+                                [](std::shared_ptr<const IMesh> p_mesh,
+                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                     bc_descriptor_list,
+                                   FunctionSymbolId rho_id, FunctionSymbolId u_id, FunctionSymbolId p_id) -> void {
+                                  switch (p_mesh->dimension()) {
+                                  case 1: {
+                                    GlaceScheme<1>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                    break;
+                                  }
+                                  case 2: {
+                                    GlaceScheme<2>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                    break;
+                                  }
+                                  case 3: {
+                                    GlaceScheme<3>{*p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                    break;
+                                  }
+                                  default: {
+                                    throw UnexpectedError("invalid mesh dimension");
+                                  }
+                                  }
                                 }
-                              }
 
-                            }));
+                              }));
 }
diff --git a/src/language/utils/PugsFunctionAdapter.hpp b/src/language/utils/PugsFunctionAdapter.hpp
index 33e844ec886c08d16e049e71a38c93476a7dd697..5e9a6cf763bb892d2b2eedf1b1c91c26c8d64f16 100644
--- a/src/language/utils/PugsFunctionAdapter.hpp
+++ b/src/language/utils/PugsFunctionAdapter.hpp
@@ -198,30 +198,89 @@ class PugsFunctionAdapter<OutputType(InputType...)>
   [[nodiscard]] PUGS_INLINE static std::function<OutputType(DataVariant&& result)>
   getResultConverter(const ASTNodeDataType& data_type)
   {
-    switch (data_type) {
-    case ASTNodeDataType::list_t: {
-      return [](DataVariant&& result) -> OutputType {
-        AggregateDataVariant& v = std::get<AggregateDataVariant>(result);
-        OutputType x;
-        for (size_t i = 0; i < x.dimension(); ++i) {
-          x[i] = std::get<double>(v[i]);
+    if constexpr (is_tiny_vector_v<OutputType>) {
+      switch (data_type) {
+      case ASTNodeDataType::list_t: {
+        return [](DataVariant&& result) -> OutputType {
+          AggregateDataVariant& v = std::get<AggregateDataVariant>(result);
+          OutputType x;
+
+          for (size_t i = 0; i < x.dimension(); ++i) {
+            std::visit(
+              [&](auto&& vi) {
+                using Vi_T = std::decay_t<decltype(vi)>;
+                if constexpr (std::is_arithmetic_v<Vi_T>) {
+                  x[i] = vi;
+                } else {
+                  throw UnexpectedError("expecting arithmetic value");
+                }
+              },
+              v[i]);
+          }
+          return x;
+        };
+      }
+      case ASTNodeDataType::vector_t: {
+        return [](DataVariant&& result) -> OutputType { return std::get<OutputType>(result); };
+      }
+      case ASTNodeDataType::bool_t: {
+        if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
+          return
+            [](DataVariant&& result) -> OutputType { return OutputType{static_cast<double>(std::get<bool>(result))}; };
+        } else {
+          throw UnexpectedError("unexpected data_type");
         }
-        return x;
-      };
-    }
-    case ASTNodeDataType::vector_t: {
-      return [](DataVariant&& result) -> OutputType { return std::get<OutputType>(result); };
-    }
-    case ASTNodeDataType::double_t: {
-      if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
-        return [](DataVariant&& result) -> OutputType { return OutputType{std::get<double>(result)}; };
-      } else {
+      }
+      case ASTNodeDataType::unsigned_int_t: {
+        if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
+          return [](DataVariant&& result) -> OutputType {
+            return OutputType(static_cast<double>(std::get<uint64_t>(result)));
+          };
+        } else {
+          throw UnexpectedError("unexpected data_type");
+        }
+      }
+      case ASTNodeDataType::int_t: {
+        if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
+          return [](DataVariant&& result) -> OutputType {
+            return OutputType{static_cast<double>(std::get<int64_t>(result))};
+          };
+        } else {
+          // If this point is reached must be a 0 vector
+          return [](DataVariant &&) -> OutputType { return OutputType{ZeroType{}}; };
+        }
+      }
+      case ASTNodeDataType::double_t: {
+        if constexpr (std::is_same_v<OutputType, TinyVector<1>>) {
+          return [](DataVariant&& result) -> OutputType { return OutputType{std::get<double>(result)}; };
+        } else {
+          throw UnexpectedError("unexpected data_type");
+        }
+      }
+      default: {
         throw UnexpectedError("unexpected data_type");
       }
-    }
-    default: {
-      throw UnexpectedError("unexpected data_type");
-    }
+      }
+    } else if constexpr (std::is_arithmetic_v<OutputType>) {
+      switch (data_type) {
+      case ASTNodeDataType::bool_t: {
+        return [](DataVariant&& result) -> OutputType { return std::get<bool>(result); };
+      }
+      case ASTNodeDataType::unsigned_int_t: {
+        return [](DataVariant&& result) -> OutputType { return std::get<uint64_t>(result); };
+      }
+      case ASTNodeDataType::int_t: {
+        return [](DataVariant&& result) -> OutputType { return std::get<int64_t>(result); };
+      }
+      case ASTNodeDataType::double_t: {
+        return [](DataVariant&& result) -> OutputType { return std::get<double>(result); };
+      }
+      default: {
+        throw UnexpectedError("unexpected data_type");
+      }
+      }
+    } else {
+      static_assert(std::is_arithmetic_v<OutputType>, "unexpected output type");
     }
   }
 };