diff --git a/euler/src/eval.rs b/euler/src/eval.rs
new file mode 100644
index 0000000..46bc5a3
--- /dev/null
+++ b/euler/src/eval.rs
@@ -0,0 +1,97 @@
+//! Traits for evaluating initial conditions, exact solutions, or boundary conditions
+
+use super::Float;
+use super::GAMMA;
+use ndarray::{azip, ArrayView, ArrayViewMut, Dimension};
+
+pub trait Evaluator<D: Dimension>: Send + Sync {
+    fn evaluate(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayViewMut<Float, D>,
+        rhou: ArrayViewMut<Float, D>,
+        rhov: ArrayViewMut<Float, D>,
+        e: ArrayViewMut<Float, D>,
+    );
+}
+
+/// Necessary to avoid specialisation of `Evaluator` trait for items
+/// that could implement both this and `Evaluator`
+#[derive(Clone, Debug)]
+pub struct EvaluatorPressureWrapper<'a, D: Dimension, E: EvaluatorPressure<D>>(
+    &'a E,
+    std::marker::PhantomData<D>,
+);
+
+impl<'a, D: Dimension, E: EvaluatorPressure<D>> EvaluatorPressureWrapper<'a, D, E> {
+    pub fn new(e: &'a E) -> Self {
+        Self(e, std::marker::PhantomData)
+    }
+}
+
+pub trait EvaluatorPressure<D: Dimension>: Send + Sync {
+    fn rho(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        out: ArrayViewMut<Float, D>,
+    );
+    fn u(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayView<Float, D>,
+        out: ArrayViewMut<Float, D>,
+    );
+    fn v(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayView<Float, D>,
+        out: ArrayViewMut<Float, D>,
+    );
+    fn p(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayView<Float, D>,
+        u: ArrayView<Float, D>,
+        v: ArrayView<Float, D>,
+        out: ArrayViewMut<Float, D>,
+    );
+}
+
+impl<'a, D: Dimension, BP: EvaluatorPressure<D>> Evaluator<D>
+    for EvaluatorPressureWrapper<'a, D, BP>
+{
+    fn evaluate(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        mut rho: ArrayViewMut<Float, D>,
+        mut rhou: ArrayViewMut<Float, D>,
+        mut rhov: ArrayViewMut<Float, D>,
+        mut e: ArrayViewMut<Float, D>,
+    ) {
+        let eva = &self.0;
+        eva.rho(t, x.view(), y.view(), rho.view_mut());
+        eva.u(t, x.view(), y.view(), rho.view(), rhou.view_mut());
+        eva.v(t, x.view(), y.view(), rho.view(), rhov.view_mut());
+        eva.p(t, x, y, rho.view(), rhou.view(), rhov.view(), e.view_mut());
+
+        azip!((rho in &rho, u in &rhou, v in &rhov, e in &mut e) {
+            let p = *e;
+            *e = p / (GAMMA - 1.0) + rho * (u*u + v*v) / 2.0;
+
+        });
+        azip!((rho in &rho, rhou in &mut rhou) *rhou /= rho);
+        azip!((rho in &rho, rhov in &mut rhov) *rhov /= rho);
+    }
+}
diff --git a/euler/src/lib.rs b/euler/src/lib.rs
index 54d83c7..b7801ab 100644
--- a/euler/src/lib.rs
+++ b/euler/src/lib.rs
@@ -6,6 +6,8 @@ use sbp::operators::{InterpolationOperator, SbpOperator2d, UpwindOperator2d};
 use sbp::utils::Direction;
 use sbp::Float;
 
+pub mod eval;
+
 mod vortex;
 pub use vortex::{vortex, VortexParameters, Vortice};
 
@@ -691,11 +693,25 @@ pub enum BoundaryCharacteristic {
     This,
     Grid(usize),
     Vortex(VortexParameters),
+    Eval(std::sync::Arc<dyn eval::Evaluator<ndarray::Ix1>>),
     // Vortices(Vec<VortexParameters>),
     Interpolate(usize, Box<dyn InterpolationOperator>),
     MultiGrid(Vec<(usize, usize, usize)>),
 }
 
+impl std::fmt::Debug for BoundaryCharacteristic {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            Self::This => write!(f, "This"),
+            Self::Grid(g) => write!(f, "Grid({})", g),
+            Self::Vortex(vp) => write!(f, "{:?}", vp),
+            Self::Eval(_) => write!(f, "Eval"),
+            Self::Interpolate(_, _) => write!(f, "Interpolate"),
+            Self::MultiGrid(m) => write!(f, "Multigrid: {:?}", m),
+        }
+    }
+}
+
 pub type BoundaryTerms<'a> = Direction<ArrayView2<'a, Float>>;
 pub type BoundaryCharacteristics = Direction<BoundaryCharacteristic>;
 
@@ -708,6 +724,7 @@ fn boundary_extractor<'a>(
         north: match &bc.north {
             BoundaryCharacteristic::This => field.south(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
+            BoundaryCharacteristic::Eval(_) => todo!(),
             BoundaryCharacteristic::Grid(_)
             | BoundaryCharacteristic::Interpolate(_, _)
             | BoundaryCharacteristic::MultiGrid(_) => panic!("Only working on self grid"),
@@ -715,6 +732,7 @@ fn boundary_extractor<'a>(
         south: match &bc.south {
             BoundaryCharacteristic::This => field.north(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
+            BoundaryCharacteristic::Eval(_) => todo!(),
             BoundaryCharacteristic::Grid(_)
             | BoundaryCharacteristic::Interpolate(_, _)
             | BoundaryCharacteristic::MultiGrid(_) => panic!("Only working on self grid"),
@@ -722,6 +740,7 @@ fn boundary_extractor<'a>(
         west: match &bc.west {
             BoundaryCharacteristic::This => field.east(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
+            BoundaryCharacteristic::Eval(_) => todo!(),
             BoundaryCharacteristic::Grid(_)
             | BoundaryCharacteristic::Interpolate(_, _)
             | BoundaryCharacteristic::MultiGrid(_) => panic!("Only working on self grid"),
@@ -729,6 +748,7 @@ fn boundary_extractor<'a>(
         east: match &bc.east {
             BoundaryCharacteristic::This => field.west(),
             BoundaryCharacteristic::Vortex(_params) => todo!(),
+            BoundaryCharacteristic::Eval(_) => todo!(),
             BoundaryCharacteristic::Grid(_)
             | BoundaryCharacteristic::Interpolate(_, _)
             | BoundaryCharacteristic::MultiGrid(_) => panic!("Only working on self grid"),
@@ -748,11 +768,6 @@ fn boundary_extract<'a>(
     match bc {
         BoundaryCharacteristic::This => seldir(field),
         BoundaryCharacteristic::Grid(g) => seldir(&fields[*g]),
-        BoundaryCharacteristic::Vortex(v) => {
-            let field = eb.unwrap();
-            vortexify(field.view_mut(), grid, v, time);
-            field.view()
-        }
         BoundaryCharacteristic::Interpolate(g, operator) => {
             let to = eb.unwrap();
             let fine2coarse = field.nx() < fields[*g].nx();
@@ -785,6 +800,24 @@ fn boundary_extract<'a>(
             }
             to.view()
         }
+        BoundaryCharacteristic::Vortex(v) => {
+            let field = eb.unwrap();
+            vortexify(field.view_mut(), grid, v, time);
+            field.view()
+        }
+        BoundaryCharacteristic::Eval(expr) => {
+            let field = eb.unwrap();
+            let (x, y) = grid;
+            let mut fiter = field.outer_iter_mut();
+            let (rho, rhou, rhov, e) = (
+                fiter.next().unwrap(),
+                fiter.next().unwrap(),
+                fiter.next().unwrap(),
+                fiter.next().unwrap(),
+            );
+            expr.evaluate(time, x, y, rho, rhou, rhov, e);
+            field.view()
+        }
     }
 }
 
@@ -864,6 +897,7 @@ impl BoundaryStorage {
         Self {
             north: match bt.north() {
                 BoundaryCharacteristic::Vortex(_)
+                | BoundaryCharacteristic::Eval(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
                     Some(ndarray::Array2::zeros((4, grid.nx())))
@@ -872,6 +906,7 @@ impl BoundaryStorage {
             },
             south: match bt.south() {
                 BoundaryCharacteristic::Vortex(_)
+                | BoundaryCharacteristic::Eval(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
                     Some(ndarray::Array2::zeros((4, grid.nx())))
@@ -880,6 +915,7 @@ impl BoundaryStorage {
             },
             east: match bt.east() {
                 BoundaryCharacteristic::Vortex(_)
+                | BoundaryCharacteristic::Eval(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
                     Some(ndarray::Array2::zeros((4, grid.ny())))
@@ -888,6 +924,7 @@ impl BoundaryStorage {
             },
             west: match bt.west() {
                 BoundaryCharacteristic::Vortex(_)
+                | BoundaryCharacteristic::Eval(_)
                 | BoundaryCharacteristic::Interpolate(_, _)
                 | BoundaryCharacteristic::MultiGrid(_) => {
                     Some(ndarray::Array2::zeros((4, grid.ny())))
diff --git a/multigrid/Cargo.toml b/multigrid/Cargo.toml
index 0370fd0..69dcd5e 100644
--- a/multigrid/Cargo.toml
+++ b/multigrid/Cargo.toml
@@ -17,3 +17,4 @@ serde = { version = "1.0.115", features = ["derive"] }
 json5 = "0.3.0"
 indexmap = { version = "1.5.2", features = ["serde-1"] }
 argh = "0.1.4"
+evalexpr = "6.2.0"
diff --git a/multigrid/examples/double.json b/multigrid/examples/double.json
index 93775af..e5295ff 100644
--- a/multigrid/examples/double.json
+++ b/multigrid/examples/double.json
@@ -5,10 +5,6 @@
             "operators": {
                 "xi": "upwind9",
                 "eta": "upwind9h2"
-            },
-            "boundary_conditions": {
-                "east": "vortex",
-                "west": "vortex"
             }
         },
         "grid0": {
@@ -26,14 +22,18 @@
             }
         }
     },
-    "integration_time": 2.0,
-    "vortex": {
-        "vortices": [{
-            "x0": -1.0,
-            "y0": 0.0,
-            "rstar": 0.5,
-            "eps": 1.0
-        }],
-        "mach": 0.5
-    }
+
+    "initial_conditions": {
+        "vortex": {
+            "vortices": [{
+                "x0": -1.0,
+                "y0": 0.0,
+                "rstar": 0.5,
+                "eps": 1.0
+            }],
+            "mach": 0.5
+        }
+    },
+    "boundary_conditions": "initial_conditions",
+    "integration_time": 2.0
 }
diff --git a/multigrid/examples/interpolation.json b/multigrid/examples/interpolation.json
index 945cc7e..22ad697 100644
--- a/multigrid/examples/interpolation.json
+++ b/multigrid/examples/interpolation.json
@@ -4,12 +4,6 @@
             "operators": {
                 "xi": "upwind9",
                 "eta": "upwind9"
-            },
-            "boundary_conditions": {
-                "south": "vortex",
-                "north": "vortex",
-                "east": "vortex",
-                "west": "vortex"
             }
         },
         "grid0": {
@@ -28,13 +22,16 @@
         }
     },
     "integration_time": 2.0,
-    "vortex": {
-        "vortices": [{
-            "x0": -1.0,
-            "y0": 0.0,
-            "rstar": 0.5,
-            "eps": 1.0
-        }],
-        "mach": 0.5
-    }
+    "initial_conditions": {
+        "vortex": {
+            "vortices": [{
+                "x0": -1.0,
+                "y0": 0.0,
+                "rstar": 0.5,
+                "eps": 1.0
+            }],
+            "mach": 0.5
+        }
+    },
+    "boundary_conditions": "initial_conditions"
 }
diff --git a/multigrid/examples/mix_operators.json b/multigrid/examples/mix_operators.json
index 4db0e80..f6a2768 100644
--- a/multigrid/examples/mix_operators.json
+++ b/multigrid/examples/mix_operators.json
@@ -4,12 +4,6 @@
             "operators": {
                 "xi": "upwind9",
                 "eta": "upwind9"
-            },
-            "boundary_conditions": {
-                "south": "vortex",
-                "north": "vortex",
-                "east": "vortex",
-                "west": "vortex"
             }
         },
         "grid0": {
@@ -56,13 +50,16 @@
         }
     },
     "integration_time": 2.0,
-    "vortex": {
-        "vortices": [{
-            "x0": -1.0,
-            "y0": 0.0,
-            "rstar": 0.5,
-            "eps": 1.0
-        }],
-        "mach": 0.5
-    }
+    "initial_conditions": {
+        "vortex": {
+            "vortices": [{
+                "x0": -1.0,
+                "y0": 0.0,
+                "rstar": 0.5,
+                "eps": 1.0
+            }],
+            "mach": 0.5
+        }
+    },
+    "boundary_conditions": "initial_conditions"
 }
diff --git a/multigrid/examples/quad.json b/multigrid/examples/quad.json
index b9c6a73..66ca823 100644
--- a/multigrid/examples/quad.json
+++ b/multigrid/examples/quad.json
@@ -48,13 +48,15 @@
         }
     },
     "integration_time": 2.0,
-    "vortex": {
-        "vortices": [{
-            "x0": -1.0,
-            "y0": 0.0,
-            "eps": 1.0,
-            "rstar": 0.5
-        }],
-        "mach": 0.5
+    "initial_conditions": {
+        "vortex": {
+            "vortices": [{
+                "x0": -1.0,
+                "y0": 0.0,
+                "eps": 1.0,
+                "rstar": 0.5
+            }],
+            "mach": 0.5
+        }
     }
 }
diff --git a/multigrid/examples/sedecim.json b/multigrid/examples/sedecim.json
index 20ceda7..b2c85e3 100644
--- a/multigrid/examples/sedecim.json
+++ b/multigrid/examples/sedecim.json
@@ -168,13 +168,15 @@
         }
     },
     "integration_time": 2.0,
-    "vortex": {
-        "vortices": [{
-            "x0": -1.0,
-            "y0": 0.0,
-            "rstar": 0.5,
-            "eps": 1.0
-        }],
-        "mach": 0.5
+    "initial_conditions": {
+        "vortex": {
+            "vortices": [{
+                "x0": -1.0,
+                "y0": 0.0,
+                "rstar": 0.5,
+                "eps": 1.0
+            }],
+            "mach": 0.5
+        }
     }
 }
diff --git a/multigrid/examples/tsection.json b/multigrid/examples/tsection.json
index d03d5d1..ee40aa4 100644
--- a/multigrid/examples/tsection.json
+++ b/multigrid/examples/tsection.json
@@ -4,12 +4,6 @@
             "operators": {
                 "xi": "upwind9",
                 "eta": "upwind9"
-            },
-            "boundary_conditions": {
-                "south": "vortex",
-                "north": "vortex",
-                "east": "vortex",
-                "west": "vortex"
             }
         },
         "grid0": {
@@ -37,13 +31,16 @@
         }
     },
     "integration_time": 2.0,
-    "vortex": {
-        "vortices": [{
-            "x0": 0.0,
-            "y0": 0.0,
-            "rstar": 0.5,
-            "eps": 1.0
-        }],
-        "mach": 0.5
-    }
+    "initial_conditions": {
+        "vortex": {
+            "vortices": [{
+                "x0": 0.0,
+                "y0": 0.0,
+                "rstar": 0.5,
+                "eps": 1.0
+            }],
+            "mach": 0.5
+        }
+    },
+    "boundary_conditions": "initial_conditions"
 }
diff --git a/multigrid/src/eval.rs b/multigrid/src/eval.rs
new file mode 100644
index 0000000..18f2cc7
--- /dev/null
+++ b/multigrid/src/eval.rs
@@ -0,0 +1,330 @@
+use euler::GAMMA;
+use evalexpr::*;
+use ndarray::{azip, ArrayView, ArrayViewMut, Dimension};
+use sbp::Float;
+
+#[derive(Clone, Debug)]
+pub enum Evaluator {
+    Pressure(EvaluatorPressure),
+    Conservation(EvaluatorConservation),
+}
+
+#[derive(Debug, Clone)]
+pub struct EvaluatorConservation {
+    pub(crate) ctx: HashMapContext,
+    pub(crate) rho: Node,
+    pub(crate) rhou: Node,
+    pub(crate) rhov: Node,
+    pub(crate) e: Node,
+}
+
+impl<D: ndarray::Dimension> euler::eval::Evaluator<D> for Evaluator {
+    fn evaluate(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayViewMut<Float, D>,
+        rhou: ArrayViewMut<Float, D>,
+        rhov: ArrayViewMut<Float, D>,
+        e: ArrayViewMut<Float, D>,
+    ) {
+        match self {
+            Self::Conservation(c) => c.evaluate(t, x, y, rho, rhou, rhov, e),
+            Self::Pressure(p) => {
+                euler::eval::EvaluatorPressureWrapper::new(p).evaluate(t, x, y, rho, rhou, rhov, e)
+            }
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct EvaluatorPressure {
+    pub(crate) ctx: HashMapContext,
+    pub(crate) rho: Node,
+    pub(crate) u: Node,
+    pub(crate) v: Node,
+    pub(crate) p: Node,
+}
+
+struct ContextWrapper<'a> {
+    ctx: &'a HashMapContext,
+    x: Option<Value>,
+    y: Option<Value>,
+    t: Value,
+    rho: Option<Value>,
+    u: Option<Value>,
+    v: Option<Value>,
+    rhou: Option<Value>,
+    rhov: Option<Value>,
+}
+
+impl<'a> ContextWrapper<'a> {
+    fn wrap(ctx: &'a HashMapContext, t: Value) -> Self {
+        Self {
+            ctx,
+            t,
+            x: None,
+            y: None,
+            rho: None,
+            rhou: None,
+            rhov: None,
+            u: None,
+            v: None,
+        }
+    }
+}
+
+impl Context for ContextWrapper<'_> {
+    fn get_value(&self, identifier: &str) -> Option<&Value> {
+        match identifier {
+            "t" => Some(&self.t),
+            "x" => self.x.as_ref(),
+            "y" => self.y.as_ref(),
+            "rho" => self.rho.as_ref(),
+            "rhou" => self.rhou.as_ref(),
+            "rhov" => self.rhov.as_ref(),
+            "u" => self.u.as_ref(),
+            "v" => self.v.as_ref(),
+            id => self.ctx.get_value(id),
+        }
+    }
+
+    fn call_function(&self, identifier: &str, argument: &Value) -> EvalexprResult<Value> {
+        self.ctx.call_function(identifier, argument)
+    }
+}
+
+impl<D: Dimension> euler::eval::Evaluator<D> for EvaluatorConservation {
+    fn evaluate(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        mut rho: ArrayViewMut<Float, D>,
+        mut rhou: ArrayViewMut<Float, D>,
+        mut rhov: ArrayViewMut<Float, D>,
+        mut e: ArrayViewMut<Float, D>,
+    ) {
+        let mut ctx = ContextWrapper::wrap(&self.ctx, t.into());
+
+        azip!((&x in &x.view(), &y in &y, rho in &mut rho) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+
+            *rho = self.rho.eval_number_with_context(&ctx).unwrap();
+        });
+
+        azip!((&x in &x, &y in &y, &rho in &rho, rhou in &mut rhou) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+            ctx.rho = Some(rho.into());
+
+            *rhou = self.rhou.eval_number_with_context(&ctx).unwrap();
+        });
+
+        azip!((&x in &x, &y in &y, &rho in &rho, rhov in &mut rhov) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+            ctx.rho = Some(rho.into());
+
+            *rhov = self.rhov.eval_number_with_context(&ctx).unwrap();
+        });
+
+        azip!((&x in &x, &y in &y, &rho in &rho, &rhou in &rhou, &rhov in &rhov, e in &mut e) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+            ctx.rho = Some(rho.into());
+            ctx.rhou = Some(rhou.into());
+            ctx.rhov = Some(rhov.into());
+
+            *e = self.e.eval_number_with_context(&ctx).unwrap();
+        });
+    }
+}
+
+impl<D: Dimension> euler::eval::EvaluatorPressure<D> for EvaluatorPressure {
+    fn rho(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        mut rho: ArrayViewMut<Float, D>,
+    ) {
+        let mut ctx = ContextWrapper::wrap(&self.ctx, t.into());
+
+        azip!((&x in &x, &y in &y, rho in &mut rho) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+
+            *rho = self.rho.eval_number_with_context(&ctx).unwrap();
+        })
+    }
+
+    fn u(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayView<Float, D>,
+        mut u: ArrayViewMut<Float, D>,
+    ) {
+        let mut ctx = ContextWrapper::wrap(&self.ctx, t.into());
+
+        azip!((&x in &x, &y in &y, &rho in &rho, u in &mut u) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+            ctx.rho = Some(rho.into());
+
+            *u = self.u.eval_number_with_context(&ctx).unwrap();
+        })
+    }
+    fn v(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayView<Float, D>,
+        mut v: ArrayViewMut<Float, D>,
+    ) {
+        let mut ctx = ContextWrapper::wrap(&self.ctx, t.into());
+
+        azip!((&x in &x, &y in &y, &rho in &rho, v in &mut v) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+            ctx.rho = Some(rho.into());
+
+            *v = self.v.eval_number_with_context(&ctx).unwrap();
+        })
+    }
+
+    fn p(
+        &self,
+        t: Float,
+        x: ArrayView<Float, D>,
+        y: ArrayView<Float, D>,
+        rho: ArrayView<Float, D>,
+        u: ArrayView<Float, D>,
+        v: ArrayView<Float, D>,
+        mut p: ArrayViewMut<Float, D>,
+    ) {
+        let mut ctx = ContextWrapper::wrap(&self.ctx, t.into());
+
+        azip!((&x in &x, &y in &y, &rho in &rho, &u in &u, &v in &v, p in &mut p) {
+            ctx.x = Some(x.into());
+            ctx.y = Some(y.into());
+            ctx.rho = Some(rho.into());
+            ctx.u = Some(u.into());
+            ctx.v = Some(v.into());
+
+            *p = self.p.eval_number_with_context(&ctx).unwrap();
+        })
+    }
+}
+
+#[test]
+fn append_default_context() {
+    let basic_ctx = context_map! {
+        "a" => 2,
+    }
+    .unwrap();
+
+    let mut ctx = ContextWrapper::wrap(&basic_ctx, Value::from(4));
+    ctx.x = Some(3.into());
+
+    let expr = "a + x + t";
+
+    let node = build_operator_tree(expr).unwrap();
+
+    assert_eq!(node.eval_with_context(&ctx).unwrap().as_int().unwrap(), 9);
+}
+
+pub fn default_context() -> HashMapContext {
+    let mut context = math_consts_context! {}.unwrap();
+
+    context.set_value("GAMMA".into(), GAMMA.into()).unwrap();
+
+    context
+        .set_function(
+            "if".into(),
+            Function::new(|arg| {
+                let arg = arg.as_tuple()?;
+                if arg.len() != 3 {
+                    return Err(error::EvalexprError::WrongFunctionArgumentAmount {
+                        expected: 3,
+                        actual: arg.len(),
+                    });
+                }
+                let b = arg[0].as_boolean()?;
+                if b.into() {
+                    Ok(arg[1].clone())
+                } else {
+                    Ok(arg[2].clone())
+                }
+            }),
+        )
+        .unwrap();
+
+    context
+        .set_function(
+            "case".into(),
+            Function::new(|arg| {
+                let mut arg = arg.as_tuple()?;
+                if arg.len() % 2 == 0 {
+                    return Err(error::EvalexprError::WrongFunctionArgumentAmount {
+                        expected: arg.len() + 1,
+                        actual: arg.len(),
+                    });
+                }
+                //let mut arg = arg.as_slice();
+                while arg.len() > 2 {
+                    let boolean = arg.remove(0);
+                    let value = arg.remove(0);
+
+                    if boolean.as_boolean()? {
+                        return Ok(value);
+                    }
+                }
+                Ok(arg.pop().unwrap())
+            }),
+        )
+        .unwrap();
+
+    context
+        .set_function(
+            "math::atan2".into(),
+            Function::new(|arg| {
+                let arg = arg.as_tuple()?;
+                if arg.len() != 2 {
+                    return Err(error::EvalexprError::WrongFunctionArgumentAmount {
+                        expected: 2,
+                        actual: arg.len(),
+                    });
+                }
+                let s = arg[0].as_number()?;
+                let o = arg[1].as_number()?;
+                Ok(s.atan2(o).into())
+            }),
+        )
+        .unwrap();
+
+    context
+        .set_function(
+            "math::hypot".into(),
+            Function::new(|arg| {
+                let arg = arg.as_tuple()?;
+                if arg.len() != 2 {
+                    return Err(error::EvalexprError::WrongFunctionArgumentAmount {
+                        expected: 2,
+                        actual: arg.len(),
+                    });
+                }
+                let s = arg[0].as_number()?;
+                let o = arg[1].as_number()?;
+                Ok(s.hypot(o).into())
+            }),
+        )
+        .unwrap();
+
+    context
+}
diff --git a/multigrid/src/main.rs b/multigrid/src/main.rs
index a4b39ef..ae3e6b0 100644
--- a/multigrid/src/main.rs
+++ b/multigrid/src/main.rs
@@ -1,12 +1,14 @@
 use argh::FromArgs;
 use rayon::prelude::*;
 
+use euler::eval::Evaluator;
 use sbp::operators::SbpOperator2d;
 use sbp::*;
 
 mod file;
 mod parsing;
 use file::*;
+mod eval;
 
 struct System {
     fnow: Vec<euler::Field>,
@@ -283,14 +285,35 @@ fn main() {
     let parsing::RuntimeConfiguration {
         names,
         grids,
-        bc: bt,
+        grid_connections,
         op: operators,
         integration_time,
-        vortex: vortexparams,
+        initial_conditions,
+        boundary_conditions: _,
     } = config.into_runtime();
 
-    let mut sys = System::new(grids, bt, operators);
-    sys.vortex(0.0, &vortexparams);
+    let mut sys = System::new(grids, grid_connections, operators);
+    match initial_conditions {
+        /*
+        parsing::InitialConditions::File(f) => {
+            for grid in &sys.grids {
+                // Copy initial conditions from file, requires name of field
+                todo!()
+            }
+        }
+        */
+        parsing::InitialConditions::Vortex(vortexparams) => sys.vortex(0.0, &vortexparams),
+        parsing::InitialConditions::Expressions(expr) => {
+            let t = 0.0;
+            for (grid, field) in sys.grids.iter().zip(sys.fnow.iter_mut()) {
+                // Evaluate the expressions on all variables
+                let x = grid.x();
+                let y = grid.y();
+                let (rho, rhou, rhov, e) = field.components_mut();
+                (*expr).evaluate(t, x, y, rho, rhou, rhov, e);
+            }
+        }
+    }
 
     let dt = sys.max_dt();
 
@@ -348,16 +371,23 @@ fn main() {
     }
 
     output.add_timestep(ntime, &sys.fnow);
+    /*
     if opt.error {
         let time = ntime as Float * dt;
         let mut e = 0.0;
         for ((fmod, grid), op) in sys.fnow.iter().zip(&sys.grids).zip(&sys.operators) {
             let mut fvort = fmod.clone();
-            fvort.vortex(grid.x(), grid.y(), time, &vortexparams);
+            fvort.vortex(
+                grid.x(),
+                grid.y(),
+                time,
+                &vortexparams.as_ref().unwrap().clone(),
+            );
             e += fmod.h2_err(&fvort, &**op);
         }
         outinfo.error = Some(e);
     }
+    */
 
     if opt.output_json {
         println!("{}", json5::to_string(&outinfo).unwrap());
diff --git a/multigrid/src/parsing.rs b/multigrid/src/parsing.rs
index 241592b..a66ae22 100644
--- a/multigrid/src/parsing.rs
+++ b/multigrid/src/parsing.rs
@@ -1,9 +1,13 @@
+use std::convert::{TryFrom, TryInto};
+
 use sbp::operators::SbpOperator2d;
 use sbp::utils::h2linspace;
 use sbp::Float;
 
 use serde::{Deserialize, Serialize};
 
+use crate::eval;
+
 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
 #[serde(rename_all = "lowercase")]
 pub enum Operator {
@@ -35,6 +39,10 @@ pub struct Linspace {
 pub enum GridLike {
     Linspace(Linspace),
     Array(ArrayForm),
+    /*
+    #[serde(rename = "initial_conditions")]
+    InitialConditions,
+    */
 }
 
 impl From<GridLike> for ArrayForm {
@@ -46,6 +54,7 @@ impl From<GridLike> for ArrayForm {
                 ndarray::Array::linspace(lin.start, lin.end, lin.steps)
             }),
             GridLike::Array(arr) => arr,
+            // GridLike::InitialConditions => Self::Unknown,
         }
     }
 }
@@ -137,25 +146,185 @@ pub struct GridConfig {
 type Grids = indexmap::IndexMap<String, GridConfig>;
 
 #[derive(Clone, Debug, Serialize, Deserialize)]
+/// Will be evaluated by evalexpr
+pub struct ExpressionsConservation {
+    pub globals: Option<String>,
+    pub rho: String,
+    pub rhou: String,
+    pub rhov: String,
+    pub e: String,
+}
+
+#[derive(Clone, Debug, Serialize, Deserialize)]
+/// Will be evaluated by evalexpr
+pub struct ExpressionsPressure {
+    pub globals: Option<String>,
+    pub rho: String,
+    pub u: String,
+    pub v: String,
+    pub p: String,
+}
+
+#[derive(Clone, Debug, Serialize, Deserialize)]
+#[serde(rename_all = "lowercase")]
+#[serde(untagged)]
+pub enum Expressions {
+    Conservation(ExpressionsConservation),
+    Pressure(ExpressionsPressure),
+}
+
+#[derive(Clone, Debug, Serialize, Deserialize)]
+#[serde(rename_all = "lowercase")]
+pub enum InputInitialConditions {
+    Vortex(euler::VortexParameters),
+    // File(String),
+    Expressions(Expressions),
+}
+
+#[derive(Clone, Debug)]
+pub enum InitialConditions {
+    Vortex(euler::VortexParameters),
+    // File(hdf5::File),
+    Expressions(std::sync::Arc<eval::Evaluator>),
+}
+
+impl TryFrom<Expressions> for eval::Evaluator {
+    type Error = ();
+    fn try_from(expr: Expressions) -> Result<Self, Self::Error> {
+        let mut context = eval::default_context();
+        match expr {
+            Expressions::Pressure(ExpressionsPressure {
+                globals,
+                rho,
+                u,
+                v,
+                p,
+            }) => {
+                if let Some(globals) = &globals {
+                    evalexpr::eval_with_context_mut(globals, &mut context).unwrap();
+                }
+                let [rho, u, v, p] = [
+                    evalexpr::build_operator_tree(&rho).unwrap(),
+                    evalexpr::build_operator_tree(&u).unwrap(),
+                    evalexpr::build_operator_tree(&v).unwrap(),
+                    evalexpr::build_operator_tree(&p).unwrap(),
+                ];
+                Ok(eval::Evaluator::Pressure(eval::EvaluatorPressure {
+                    ctx: context,
+                    rho,
+                    u,
+                    v,
+                    p,
+                }))
+            }
+            Expressions::Conservation(ExpressionsConservation {
+                globals,
+                rho,
+                rhou,
+                rhov,
+                e,
+            }) => {
+                if let Some(globals) = &globals {
+                    evalexpr::eval_with_context_mut(globals, &mut context).unwrap();
+                }
+                let [rho, rhou, rhov, e] = [
+                    evalexpr::build_operator_tree(&rho).unwrap(),
+                    evalexpr::build_operator_tree(&rhou).unwrap(),
+                    evalexpr::build_operator_tree(&rhov).unwrap(),
+                    evalexpr::build_operator_tree(&e).unwrap(),
+                ];
+                Ok(eval::Evaluator::Conservation(eval::EvaluatorConservation {
+                    ctx: context,
+                    rho,
+                    rhou,
+                    rhov,
+                    e,
+                }))
+            }
+        }
+    }
+}
+
+impl TryFrom<InputInitialConditions> for InitialConditions {
+    type Error = ();
+    fn try_from(v: InputInitialConditions) -> Result<Self, Self::Error> {
+        Ok(match v {
+            InputInitialConditions::Vortex(v) => Self::Vortex(v),
+            // InputInitialConditions::File(file) => Self::File(hdf5::File::open(file).unwrap()),
+            InputInitialConditions::Expressions(expr) => {
+                Self::Expressions(std::sync::Arc::new(expr.try_into()?))
+            }
+        })
+    }
+}
+
+#[derive(Clone, Debug, Serialize, Deserialize)]
+#[serde(rename_all = "lowercase")]
+pub enum InputBoundaryConditions {
+    /// Initial conditions also contain the bc
+    #[serde(rename = "initial_conditions")]
+    InputInitialConditions,
+    Vortex(euler::VortexParameters),
+    Expressions(Expressions),
+    #[serde(rename = "not_needed")]
+    NotNeeded,
+}
+
+impl Default for InputBoundaryConditions {
+    fn default() -> Self {
+        Self::NotNeeded
+    }
+}
+
+#[derive(Clone, Debug)]
+pub enum BoundaryConditions {
+    Vortex(euler::VortexParameters),
+    Expressions(std::sync::Arc<eval::Evaluator>),
+    NotNeeded,
+}
+
+#[derive(Clone, Debug, Serialize, Deserialize)]
+/// Input configuration (json)
 pub struct Configuration {
     pub grids: Grids,
     pub integration_time: Float,
-    pub vortex: euler::VortexParameters,
+    pub initial_conditions: InputInitialConditions,
+    #[serde(default)]
+    pub boundary_conditions: InputBoundaryConditions,
 }
 
 pub struct RuntimeConfiguration {
     pub names: Vec<String>,
     pub grids: Vec<sbp::grid::Grid>,
-    pub bc: Vec<euler::BoundaryCharacteristics>,
+    pub grid_connections: Vec<euler::BoundaryCharacteristics>,
     pub op: Vec<Box<dyn SbpOperator2d>>,
     pub integration_time: Float,
-    pub vortex: euler::VortexParameters,
+    pub initial_conditions: InitialConditions,
+    pub boundary_conditions: BoundaryConditions,
 }
 
 impl Configuration {
     pub fn into_runtime(mut self) -> RuntimeConfiguration {
         let default = self.grids.shift_remove("default").unwrap_or_default();
         let names = self.grids.keys().cloned().collect();
+
+        let initial_conditions: InitialConditions =
+            self.initial_conditions.clone().try_into().unwrap();
+
+        let boundary_conditions = match &self.boundary_conditions {
+            InputBoundaryConditions::NotNeeded => BoundaryConditions::NotNeeded,
+            InputBoundaryConditions::Vortex(vp) => BoundaryConditions::Vortex(vp.clone()),
+            InputBoundaryConditions::Expressions(expr) => BoundaryConditions::Expressions(
+                std::sync::Arc::new(expr.clone().try_into().unwrap()),
+            ),
+            InputBoundaryConditions::InputInitialConditions => match &initial_conditions {
+                InitialConditions::Vortex(vp) => BoundaryConditions::Vortex(vp.clone()),
+                InitialConditions::Expressions(expr) => {
+                    BoundaryConditions::Expressions(expr.clone())
+                } // _ => panic!("Boundary conditions were set to initial conditions, although initial conditions are not available",),
+            },
+        };
+
         let grids = self
             .grids
             .iter()
@@ -198,7 +367,23 @@ impl Configuration {
                     (ArrayForm::Array2(x), ArrayForm::Array2(y)) => {
                         assert_eq!(x.shape(), y.shape());
                         (x, y)
-                    }
+                    } /*
+                      (ArrayForm::Unknown, ArrayForm::Unknown) => {
+                          if let InitialConditions::File(file) = &initial_conditions {
+                              let g = file.group(name).unwrap();
+                              let x = g.dataset("x").unwrap().read_2d::<Float>().unwrap();
+                              let y = g.dataset("y").unwrap().read_2d::<Float>().unwrap();
+                              assert_eq!(x.shape(), y.shape());
+                              (x, y)
+                          } else {
+                              panic!(
+                                  "Grid {} requires a valid file for setting initial size",
+                                  name
+                              );
+                          }
+                      }
+                      _ => todo!(),
+                      */
                 };
                 sbp::grid::Grid::new(x, y).unwrap()
             })
@@ -237,11 +422,11 @@ impl Configuration {
                 Box::new((matcher(eta), matcher(xi))) as Box<dyn SbpOperator2d>
             })
             .collect();
-        let bc = self
+        let grid_connections = self
             .grids
             .iter()
             .enumerate()
-            .map(|(i, (_name, g))| {
+            .map(|(i, (name, g))| {
                 let default_bc = default.boundary_conditions.clone().unwrap_or_default();
                 g.boundary_conditions
                     .clone()
@@ -249,10 +434,26 @@ impl Configuration {
                     .zip(default_bc)
                     .map(|(bc, fallback)| bc.or(fallback))
                     .map(|bc| match bc {
-                        None | Some(BoundaryType::Vortex) => {
-                            euler::BoundaryCharacteristic::Vortex(self.vortex.clone())
-                        }
+                        None => match &boundary_conditions {
+                            BoundaryConditions::Vortex(vortex) => {
+                                euler::BoundaryCharacteristic::Vortex(vortex.clone())
+                            }
+                            BoundaryConditions::Expressions(expr) => {
+                                euler::BoundaryCharacteristic::Eval(expr.clone() )
+                            }
+                            _ => panic!(
+                                "Boundary conditions are not available, but needed for grid {}",
+                                name
+                            ),
+                        },
                         Some(BoundaryType::This) => euler::BoundaryCharacteristic::Grid(i),
+                        Some(BoundaryType::Vortex) => euler::BoundaryCharacteristic::Vortex(
+                            if let BoundaryConditions::Vortex(vortex) = &boundary_conditions {
+                                vortex.clone()
+                            } else {
+                                panic!("Wanted vortex boundary conditions not found, needed for grid {}", name)
+                            },
+                        ),
                         Some(BoundaryType::Neighbour(name)) => {
                             let j = self.grids.get_index_of(&name).unwrap();
                             euler::BoundaryCharacteristic::Grid(j)
@@ -282,10 +483,11 @@ impl Configuration {
         RuntimeConfiguration {
             names,
             grids,
-            bc,
+            grid_connections,
             op,
             integration_time: self.integration_time,
-            vortex: self.vortex,
+            initial_conditions,
+            boundary_conditions,
         }
     }
 }
@@ -298,6 +500,11 @@ pub enum ArrayForm {
     Array1(ndarray::Array1<Float>),
     /// The usize is the inner dimension (nx)
     Array2(ndarray::Array2<Float>),
+    /*
+    /// A still unknown array, will be filled out by later
+    /// pass when initial_conditions file is known
+    Unknown,
+    */
 }
 
 impl From<ndarray::Array1<Float>> for ArrayForm {