Below is given the program run on each one of the several thousand simulated processors. It implements the momentum-conserving model.
;----------------------------------------------------------------------------
; transfer.scm
; Conservative exchange process for simulating physical phenomena
; Code for an individual processor.
; See HLSIM documentation
;----------------------------------------------------------------------------
(declare (usual-integrations))
;; The function exchange-func defines the interaction that takes place
;; between a pair of processors. For processors i and j, it takes
;; q_i p_i q_j p_j and returns a list of their new values.
;; The total amount of the conserved quantity in the pair remains the
;; same.
;; Diffusion:
;; This function implements diffusion of a substance represented by q,
;; and conserves q. The p variables are ignored.
(define (diffusion-exchange-func q_i p_i q_j p_j)
(let ((delta-q (* k (- q_j q_i))))
(list (+ q_i delta-q) 0.0 (- q_j delta-q) 0.0)
)
)
;; Wave:
;; This function conserves momentum.
(define (momentum-conserving-exchange-func q_i p_i q_j p_j)
(let* ((delta-p (* k (- q_j q_i) delta-t))
(p_i-new (+ p_i delta-p))
(p_j-new (- p_j delta-p)))
(list (+ q_i (* p_i-new delta-t)) p_i-new
(+ q_j (* p_j-new delta-t)) p_j-new)
)
)
(define exchange-func momentum-conserving-exchange-func)
;; the "spring constant"
(define k 1.0); 50.0)
(define delta-t 0.01)
;; This is the fraction that the time between exchanges varies by.
(define variation ;; 0.1 good for 1d, 0.05 for 2d
(if one-d? 0.1 0.05))
;----------------------------------------------------------------------------
;; period defines the average time between exchanges
(define period 1000)
;; To get a display for 1d where q is plotted on a vertical axis, set this.
;; Otherwise, a range of colors (blue = -1, white = 0, black = 1) will be used
;; to indicate q.
(define show-heights? one-d?)
;; Parameters derived from the above
(define time-delta (* variation period))
(define base-period (- period (* period (/ variation 2))))
(define (timing-func) (+ base-period (random time-delta)))
(define show-self
(if show-heights? show-height shade2))
(define show-neighbor
(if show-heights? show-neighbor-height shade-neighbor2))
;----------------------------------------------------------------------------
; The exchange process
;; does the exchange between this processor and another, and displays the
;; new values.
(define (do-exchange other)
(let* ((current-p (my-p)) (current-neighbor-p (neighbor-p other))
(current-q (my-q)) (current-neighbor-q (neighbor-q other))
(result (momentum-conserving-exchange-func current-q current-p
current-neighbor-q current-neighbor-p))
)
(set-my-q! (enforce-boundary-conditions (first result)))
(set-my-p! (second result))
(set-neighbor-q! other
(enforce-neighbor-boundary-conditions other (third result)))
(set-neighbor-p! other (fourth result))
(show-self (my-q) current-q)
(show-neighbor other (neighbor-q other) current-neighbor-q)
)
)
;; a process that periodically exchanges with a given neighbor, used below
(define (wait-and-do-exchange neighbor)
(define (loop)
(pause (timing-func))
(externally (lambda () (do-exchange neighbor)))
(loop)
)
(loop)
)
;; This sets up a thread for each neighbor. Each thread pauses and then
;; does the exchange with that neighbor. The start times of the threads
;; are staggered. The effect is to exchange with neighbors one at a time,
;; but the order changes over time.
(define (activity)
(define (start lst)
(if (null? lst)
'()
(begin
(make&start-thread (lambda () (pause (random period))
(wait-and-do-exchange (car lst))))
(start (cdr lst)))
)
)
(start neighbor-list)
)
;; This exchanges with neighbors one at a time
(define (activity2)
(pause (timing-func))
(externally (lambda () (do-exchange (random-neighbor))))
(activity2)
)
;----------------------------------------------------------------------------
; Boundary conditions
;; For fixed boundaries: returns #t if particle n is on the boundary
(define (on-boundary? n)
(or (< (neighbor-x n) 0.005) (> (neighbor-x n) 0.995)
(< (neighbor-y n) 0.005) (> (neighbor-y n) 0.995))
;; for a parabola:
;; (> (neighbor-x n) (- 1 (expt (* 1.3 (- .5 (neighbor-y n))) 2)))
)
(define self-on-boundary? (on-boundary? (processor-number)))
;; The value to fix the boundary at
(define boundary-value 0.0)
(define (enforce-boundary-conditions x)
(if self-on-boundary?
boundary-value
x)
)
(define (enforce-neighbor-boundary-conditions n x)
(if (on-boundary? n)
boundary-value
x)
)
;----------------------------------------------------------------------------
; Supporting functions
(define my-id (processor-number)) ;(random 1000000))
(define neighbors
(vector-ref (simulation.neighbours (the-simulation))
my-id)
)
(define neighbor-list (vector->list neighbors))
(define num-neighbors (vector-length neighbors))
(define (pause t) (wait (make-timeout-event t)))
(define (my-q) (vector-ref (simulation.get (the-simulation) 'q) my-id))
(define (my-p)
(vector-ref (simulation.get (the-simulation) 'p) my-id))
(define (neighbor-q n)
(vector-ref (simulation.get (the-simulation) 'q) n))
(define (neighbor-p n)
(vector-ref (simulation.get (the-simulation) 'p) n))
(define (set-my-q! x)
(vector-set! (simulation.get (the-simulation) 'q) my-id x))
(define (set-my-p! x)
(vector-set! (simulation.get (the-simulation) 'p) my-id x))
(define (set-neighbor-q! n x)
(vector-set! (simulation.get (the-simulation) 'q) n x))
(define (set-neighbor-p! n x)
(vector-set! (simulation.get (the-simulation) 'p) n x))
;; Get coordinates - used only in setting the initial conditions
(define (my-x) (vector-ref (simulation.x (the-simulation)) my-id))
(define (my-y) (vector-ref (simulation.y (the-simulation)) my-id))
(define (neighbor-x n) (vector-ref (simulation.x (the-simulation)) n))
(define (neighbor-y n) (vector-ref (simulation.y (the-simulation)) n))
;----------------------------------------------------------------------------
; Display code
;; Functions to display sites. There are two versions of each: one for
;; the current site and one for the neighbor with which it exchanged.
(define (shade d)
(%color-me (the-simulation) my-id
(shade-color (shade-color-number d)))
)
(define (shade-neighbor n d)
(%color-me (the-simulation) n (shade-color (shade-color-number d)))
)
(define (show-height d old)
(%erase-at (the-simulation) (my-x) (* 0.5 (+ old 1)))
(%show-at-color (the-simulation) (my-x) (* 0.5 (+ d 1)) my-color)
)
(define (show-neighbor-height n d old)
(%erase-at (the-simulation) (vector-ref (simulation.x (the-simulation)) n)
(* 0.5 (+ old 1)))
(%show-at (the-simulation) (vector-ref (simulation.x (the-simulation)) n)
(* 0.5 (+ d 1)))
)
(define my-color (if self-on-boundary? "green" "blue"))
;---------------------------------------------------------------------------
; Initial conditions
(define 2pi (* 8 (atan 1 1)))
(define pi (/ 2pi 2))
;; A sine-shaped perturbation
(define (sine-wave-initial-conditions)
(set-my-q!
(let* ((x (my-x)) (y (my-y))
(r (sqrt (+ (expt (- x .5) 2) (expt (- y .5) 2)))))
(if (< r .3)
(* (if one-d? 1.0 2.0) (+ .5 (* .5 (cos (* r (/ 2pi .6))))))
0.0)
)
)
(set-my-p! 0.0)
)
;; Standing wave (1d)
(define (standing-wave-initial-conditions)
(set-my-q!
(let* ((x (my-x)) (y (my-y))
(r (sqrt (+ (expt (- x .5) 2) (expt (- y .5) 2)))))
(* 1.0 (sin (* (/ x (/ 100 101)) pi)))
)
)
(set-my-p! 0.0)
)
(define initial-conditions standing-wave-initial-conditions)
;;=====================================================================
; Start the processor
(initial-conditions)
(activity)
(wait (make-event))