Workbook

# Proof of Concept TL;DR the results table on the right are live results of running each code block in this document, and if you press "View" on the code block at the bottom that evaluates to a `cfd/fvm.ScenarioDefinition` then it will run it live in the browser. This is a sample CFD workbook I've written up to showcase the core idea behind this project: specifying CFD scenarios in a literate lisp and web-based environment. This was done as a "walking skeleton" - ie each component (the Lisp interpreter, the CFD library and the web text editor) is as bare bones as it can be to get to a point where the core loop of edit workbook, evaluate output, view CFD, is achievable to function as a proof of concept. All that to say, it's pretty basic but better than just using your imagination. ## Brief intro to what you're looking at and some Lisp This is a "workbook", which is just plaintext that can be interpreted by an interpreter. The core premise is by default, you can just write markdown prose (which is like a very basic HTML/LaTeX) and the intepreter ignores it, but if it sees the following (ie the j and then an opening parentheses) then it interprets it as lisp code and returns the result in the results table. As this is the first code block it's the first in the table. j(define this-is "a code block") We can use these to define symbols (ie `this-is` now evaluates to the string "a code block") or procedures (functions) with the following syntax: j(define ^2 (lambda (x) (* x x))) j(^2 5) ;; this should be 25 in the results table! A `lambda` expression takes a list of parameters as it's first argument and evaluates to the result of its second argument. The rule for evaluation is very simple: if you have a list enclosed by (parentheses) then to evaluate it, the first item is the procedure and the rest are passed to the procedure as arguments. The only exceptions are special forms, like `define` which binds a symbol to some other expression, `if` which takes a predicate (test) as its first argument and then replaces itself with the first expression after that if true or the next if false, and `lambda` which returns an anonymous procedure. As a shorthand for `(define symbol (lambda (arg) expr))` you can write bound procedures like this: j(define (factorial n) (if (= n 1) 1 (* n (factorial (- n 1))))) j(factorial 10) I'll get to some CFD now but the key points are: - You can write code like a report or document - Everything evaluates to something - The language features are minimal but very powerful ## Meshing The first step of creating a fully defined scenario. Currently the meshing capabilities are limited to generating a simple uniform cartesian mesh with no internal geometry. j(define mesh (structured-mesh :nx 50 :ny 25 :width 1 :height 0.5)) Notice the keyword arguments (`kwargs`) which start with a colon. These aren't strictly a feature of Scheme, the lisp I am mostly basing this dialect on, but for a config DSL they seemed like a really nice bit of syntactic sugar. ## Fields A field (when resolved) is just a container for cell data like cellValues, cellValues0 (values from previous timestep) and faceValues. j(define dye-field (prognostic-scalar-field "dye" :initial-value 0)) Until I get a SIMPLE solver fully working, I'm using a derived vector field that takes a callback function to calculate it's uniform value every timestep. I'm going to make it sinusoidal so I'll import my maths library. j(import "maths") j(define (velocity-time-unsteadiness time) [0.3 (* 0.1 (maths/sin time))]) j(define velocity-field (derived-vector-field "velocity" velocity-time-unsteadiness)) It's not super obvious but I happen to know (as I wrote the code) that the function `derived-vector-field` "binds" to a compiled go function that takes the name of the new vector field and a function that accepts a time value and returns a vector (what the square brackets represent, notice there's no need for commas - everything is delimited by spaces) for x and y values. In the hot loop of the simulation the compiled go code will run this interpreted function to update the field's uniform value! ## Equations An `equation` is the orchestrator for assembling the system of linear equations to be solved. It takes a number of operators which themselves contain behaviour for what coefficients to add where based on the fields passed to the operator. This approach is based on the way OpenFOAM writes equations in an almost symbolic maths way. j(define density 1) ;; conserving dye concentration so density = 1 j(define diffusivity 1e-6) ;; To keep the Peclet number high j(define dye-equation (equation dye-field (ddt dye-field density) (div dye-field density velocity-field) (laplacian dye-field diffusivity) (scalar-point-src [-1 0] 0.1))) Note the point source coordinates are in [-1, 1] space rather than the actual physics space for consistency. ### Boundary conditions The equation is responsible for boundary conditions rather than the field as I want to keep the `field` data type as dumb as possible. Boundary conditions are applied with the source terms as "decorators" of the matrix and rhs, using cached flux values from flux operators corresponding to boundary connections. j(set-boundary-conditions dye-equation mesh :north-border (scalar-neumann 0) :east-border (scalar-outflow) :south-border (scalar-neumann 0) :west-border (scalar-dirichlet 0)) It's the mesh that's responsible for exposing what named boundaries need to have boundary conditions applied. If you remove one of the boundary conditions the scenario won't resolve and I don't have a solid error propogation system from compiled code back to interpreted land so for now trust that the keyword boundary conditions above are the ones required for the mesh. ## Solver This is pretty self explanatory, the solver used by a scenario. I've implemented the CG method with a basic Jacobi preconditioner as a sweeet spot for simplicity and performance on the browser, however for more advanced use cases of running many convergence cases back-to-back I'd implement a more complicated multi-threaded solver and then serialise the problem definition (ie the workbook) and send it from the browser to a server somewhere to run it at full speed. For small visualisations on the browser though CG seems pretty good. j(define solver (jacobi-preconditioned-cg :max-iterations 100 :tolerance 1e-3)) ## Scenario A `Scenario` is the main orchestrator of fields and equations, and is responsible for orchestrating an iteration step and relaying if it's reached convergence for a given time step. For a pressure-velocity solver this'll involved multiple iterations and lots of complicated bits like organising Rhie-Chow MWI for the velocity face values and caching that in an iteration etc. For my basic passive transport case I need none of that - I just assemble equations and solve and assume convergence in one step. Furthermore, a `Scenario` is responsible for the complete "resolution" of the problem - the data structures declared earlier (meshes, fields, equations) are very sparse definition types under the hood that do not have any of the long data arrays initialised yet, to avoid allocating memory while defining the problem. When it comes time to run a visualisation or run a scenario to convergence the scenario will resolve all its components and only then will the long data arrays be initialised. j(define scenario (passive-transport :mesh mesh :solver solver :fields [dye-field velocity-field] :equations [dye-equation])) If you go to the results table there should be a button to "View" the scenario definition. Please feel free to tweak and change values and see if you can get other cool behaviour - currently I don't have a way to persist changes: on refresh the page will just reset to the original workbook I've hardcoded in so go nuts. I hope it's clear that with a bit of syntax highlighting and cleverness (ie rendering the prose in a normal font and rendering headings like headings etc) this could be a really powerful way of describing fully reproducible CFD setups and more.