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

[Panicz Maciej Godek <godek.maciek@gmail.com>]

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

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 <arne_bab@web.de>:

> 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: Type: text/html, Size: 8430 bytes --]