Add Evaluator for defining IC/BC in json config

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
+}

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());

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 {