[Gsll-devel] example of gsll use

[Mirko Vukovic]

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