Liam,
here is an example of use of gsll. I use the lu matrix solvers to
calculate the radiative heat flux between three infinite parallel
plates. This function takes advantage of letm to return a *closure*
that contains several solutions to the problem.
To access the results of the calculation do the following
> (funcall *closure* args) where args is one of :mat, :per, :rhs, :q1, :fluxes
There are probably other ways of accomplishing similar things (like
values) but I like closure because it provides permanence to the
calculation until I decide to discard it.
The code is not super-clean. I wrote it a month or two back, and
there is some bit-rot in it.
You can see that I had to use the data statement to revive the data.
Otherwise, the closure was not able to access them. It is for this
reason that I asked if we could have letm not destroying by default
all of its contents, but saving some for closure use. One could label
such variables with a special symbol (like !S like foo!S, which the
macro would detect and then decide whether to save or not).
Mirko
(defun 3parallel-plate-equilibrium-flux (eps0 T0 eps1-0 eps1-2 T1 eps2 T2)
"Solve flux to middle of three parallel plates with diffuse emission
Parameters:
eps0, T0: emissivity and temperature of plate 0
eps1-0, eps1-2, T1: surface emissivities of plate 1 facing plate 1 and
2 and temperature
eps2, T2: emissivity and temperature of plate 2
Get results as follows:
q1: (send *prob* :q1): desired flux
fluxes: (send *prob* :fluxes): all fluxes
Need macro send from my-utils. Otherwise use (funcall closer args)"
(letm ((dim 5)
(mat (make-array (list dim dim)
:element-type 'double-float
:initial-contents
(list (list 1d0 (1- eps0) 0d0 0d0 0d0)
(list (1- eps1-0) 1d0 0d0 0d0 (0- eps1-0))
(list 0d0 0d0 1d0 (1- eps1-2) (0- eps1-2))
(list 0d0 0d0 (1- eps2) 1d0 0d0)
(list eps1-0 0d0 0d0 eps1-2 1d0))))
(mmat (matrix-double-float mat))
(mmat0 mmat)
(rhs (make-array dim
:element-type 'double-float
:initial-contents (list (st4 T0 eps0) 0d0 0d0
(st4 T2 eps2)
(st4 T1 (+ eps1-0 eps1-2)))))
(mrhs (vector-double-float rhs))
(per (permutation dim))
(0-vec (vector-double-float (make-array dim
:initial-element 0d0)))
(x (vector-double-float dim)))
(lu-decomp mmat per)
(lu-solve mmat per mrhs x)
(data x)
(data mmat)
(data per)
(data mrhs)
(data mmat0)
(data 0-vec)
(lambda (command)
(case command
;; (:solve (lu-solve mmat per mrhs x))
;; (:check (gemv :notrans 1d0 mmat0 x 0d0 0-vec))
(:mat mat)
(:per (data per))
(:rhs (data mrhs))
(:q1 (aref (data x) 4))
(:fluxes (data x))))))
