I wrote fluid advection code: How to make this more elegant?

[Arne Babenhauserheide <arne_bab@web.de>]

[-- Attachment #1: Type: text/plain, Size: 4745 bytes --]

Hi,

I just recreated a fluid advection exercise in Guile Scheme and I’m not
quite happy with its readability. Can you help me improve it?

My main gripe is that the math does not look instantly accessible.

The original version was in Python:

    psi[i] - c1*(psi[i+1] - psi[i-1]) + c2*(psi[i+1] - 2.0*psi[i] + psi[i-1])

My port to Scheme looks like this:

    (let ((newvalue (+ (- (psir i)
                          (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
                       (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
                                (psir (- i 1)))))))
      (array-set! psinew newvalue i))


Liebe Grüße,
Arne

Here’s the full code:

#!/bin/sh
# -*- scheme -*-
exec guile -e '(@@ (advection) main)' -s "$0" "$@"
!#

; Copyright (c) 2015 John Burkardt (original Python), 2016 Corinna
; Hoose (adaption) and 2016 Arne Babenhauserheide (pep8 + Scheme
; version).

; License: LGPL, built on the Python version from 2015 John Burkardt
; and Corinna Hoose. License LGPL.

(define-module (advection)
  #:use-module (ice-9 optargs) ; define*
  #:use-module (srfi srfi-1) ; iota
  #:use-module (ice-9 format)
  #:use-module (ice-9 popen))


(define* (fd1d-advection-lax-wendroff #:key (nx 101) (nt 1000) (c 1))
  (let* ((dx (/ 1 (- nx 1)))
         (x (iota nx 0 (/ 1 nx)))
         (dt (/ 1 nt))
         (c1 (* #e0.5 (* c (/ dt dx))))
         (c2 (* 0.5 (expt (* c (/ dt dx)) 2))))
    (format #t "CFL condition: dt (~g) ≤ (~g) dx/c\n" dt (/ dx c))
    (let ((psi (make-array 0 nx))
          (X (make-array 0 nx (+ nt 1)))
          (Y (make-array 0 nx (+ nt 1)))
          (Z (make-array 0 nx (+ nt 1))))
      (let ((psinew (let ((pn (make-array 0 nx)))
                      (let loop ((i 0))
                        (cond ((= i nx)
                               pn)
                              (else
                               (let ((xi (list-ref x i)))
                                 (when (and (<= 0.4 xi) (<= xi 0.6))
                                   (array-set! pn
                                               (* (expt (- (* 10 xi) 4) 2)
                                                  (expt (- 6 (* 10 xi)) 2))
                                               i))
                                 (loop (+ 1 i)))))))))
        (define (psir i) (array-ref psi i))
        (let loop ((j 0))
          (cond
           ((> j nt) #t) ; done
           (else 
            (let ((t (/ j nt)))
              (when (>= j 1)
                (let ((newvalue (+ (- (psir 0)
                                      (* c1 (- (psir 1)
                                               (psir (- nx 1)))))
                                   (* c2 (+ (- (psir 1)
                                               (* 2 (psir 0)))
                                            (psir (- nx 1)))))))
                  (array-set! psinew newvalue 0))
                (let loop ((i 1))
                  (when (< i (- nx 1))
                    (let ((newvalue (+ (- (psir i)
                                          (* c1 (- (psir (+ i 1)) (psir (- i 1)))))
                                       (* c2 (+ (- (psir (+ i 1)) (* 2 (psir i)))
                                                (psir (- i 1)))))))
                      (array-set! psinew newvalue i))
                    (loop (+ i 1))))
                (let ((newvalue (+ (- (psir (- nx 1))
                                      (* c1 (- (psir 0) (psir (- nx 2)))))
                                   (* c2 (+ (- (psir 0) (* 2 (psir (- nx 1))))
                                            (psir (- nx 2)))))))
                  (array-set! psinew newvalue (- nx 1))))
              (let loop ((i 0))
                (when (< i nx)
                  (array-set! psi (array-ref psinew i) i)
                  (array-set! X (list-ref x i) i j)
                  (array-set! Y t i j)
                  (array-set! Z (psir i) i j)
                  (loop (+ i 1)))))
            (loop (+ j 1)))))
      (list X Y Z)))))

(define (main args)
  (display "Hello World!\n")
  (let ((res (fd1d-advection-lax-wendroff)))
    (let ((X (list-ref res 0))
          (Y (list-ref res 1))
          (Z (list-ref res 2)))
      ; (display Z)
      (let ((port (open-output-pipe "python")))
        (format port "import pylab as pl\n")
        (format port "y = [[float(j) for j in i] for i in eval('~A'[2:].replace(' ',', '))]\n" Z)
        (format port "pl.imshow(y)\n")
        (format port "pl.colorbar()\n")
        (format port "pl.show()\n")
        (close-pipe port))))
  (newline))


-- 
Unpolitisch sein
heißt politisch sein
ohne es zu merken

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 298 bytes --]