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

[Panicz Maciej Godek]

Hi,

although I cannot be of immediate help with the topic, because I don't know anything about advection, I think that the main problem with your code is that it is imperative. While Python's stylistics handles the imperative code quite nicely, it looks rather strange in Scheme.

From my experience, the proper Scheme solution should resemble the mathematical formulation of a problem (except that it should be more descriptive), rather than a list of steps for solving it.

Also, it should be more readable to use patern matching instead of list indexing, so most likely your _expression_

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

should look like this

(map (lambda (prev this next)

          (- this
             (* c1 (- next prev))

             (* (- c2) (+ next (* -2 this) prev))))

       `(0 ,@(drop psir 1))

       psir

       `(,@(drop-right psir 1) 0))

While this may also look slightly difficult to read (and write), this isn't solely because of the _expression_'s structure, but because the factors of the _expression_ have no name, and therefore the source code doesn't explain their role (this is the problem of the Python code either, but Python doesn't prompt to fix that)

HTH

PS I think that this subject fits better to guile-user

2016-01-23 11:00 GMT+01:00 Arne Babenhauserheide <address@hidden>:

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