(in-package :clf)

;; функция-добавок к давлению
(defun sigma (phi)
  phi *sigma*)

;; производная функции sigma
(defun Dsigma (phi)
  phi 0)

;; собственные числа ячейки
;; возвращает одно из 3 значений, в зависимости от выбранного type: 'A, 'S или 'Afull
(defun eigen-values-by-cell (cell type)
  (let* ((phi (phi cell))
	 (rg (rg cell))
	 (ug (ug cell))
	 (ul (ul cell))
	 (rl *rl*)
	 (cg (/ *Cg* (sqrt *kappa*))))
      (case type

	((A flux) ;; собственные числа якобиана потоковых членов
	 (let ((cl (sqrt (/ (+ (* rg cg cg)
			       (* phi (- (sigma phi)
					 rl ul ul)))
			    (* phi rl)))))
	   (let ((x1 (- ug cg))
		 (x2 (+ ug cg))
		 (x3 (+ ul cl))
		 (x4 (- ul cl)))
	     (when (>= *debug-level* 7)
	       (print-content t cell) (newline)
	       (format t "eigens-A: ~@{~d ~}~%~%" x1 x2 x3 x4))
	     (list (abs x1) (abs x2) (abs x3) (abs x4)))))

	((S source) ;; собственные числа (одно) якобиана источниковых членов
	 (let* ((ug-ul (if (= ug ul)
			   1.d-8
			   (- ug ul)))
		(x1 (- (* (abs ug-ul)
			  (+ (/ rl
				(* phi rg))
			     (/ 1
				(- 1 phi)))))))
	   (when (>= *debug-level* 6)
	     (print-content t cell) (newline)
	     (format t "eigens-S: ~d~%~%" x1))
	   (list (abs x1))))

	((Afull) ;; собственные числа якобиана потоковых членов, объединённых с
		 ;; источниками, обусловленными сужением трубки тока
	 (labels ((frolov-equation (x)
		    (- (* phi
			  (- (* rl (expt (- x ul) 2))
			     (- (sigma phi)
				(* (- 1 phi) (Dsigma phi))))
			  (- (expt (- x ug) 2)
			     (expt cg 2)))
		       (* (expt cg 2)
			  (- 1 phi)
			  rg
			  (expt (- x ug) 2))))
		  (frolov-equation-derivative (x)
		    (- (+ (* 2 phi rl
			     (- x ul)
			     (- (expt (- x ug) 2)
				(expt cg 2)))
			  (* 2 phi
			     (- x ug)
			     (- (* rl (expt (- x ul) 2))
				(- (sigma phi)
				   (* (- 1 phi) (Dsigma phi))))))
		       (* 2 rg
			  (expt cg 2)
			  (- 1 phi)
			  (- x ug))))
		  (frolov-equation-fdf (x)
		    (values (frolov-equation x)
			    (frolov-equation-derivative x)))
		  (find-root (&optional initial (method +newton-fdfsolver+))
		    (when *debug-newton*
		      (format t "~%cell:~%")
		      (print-content t cell)
		      (newline))
		    (let* ((max-iter *newton-max-iter*)
			   (solver (make-one-dimensional-root-solver-fdf
				    method #'frolov-equation #'frolov-equation-derivative #'frolov-equation-fdf initial)))
		      (loop for iter from 0
			 for oldroot = initial then root
			 for root = (progn (iterate solver) (solution solver))
			 for olderror = 0 then error
			 for error = (- root oldroot)
			 while (if (< iter max-iter)
				   (not (root-test-delta root oldroot *newton-est-error* 0.0d0))
				   (error "Метод Ньютона превысил допустимое число итераций."))
			 do
			   (when *debug-newton*
			     (when (zerop iter)
			       (format t "NEWTON: initial-root ~d~%" oldroot))
			     (format t "NEWTON: iter ~d, root ~d, error ~d~%"
				     (1+ iter) root error))
			   (when (minusp (* olderror error))
			     (error "Нарушена монотонность метода Ньютона. Возможно, получилось отрицательное давление?"))
			 finally (return root)))))
	   (let* ((suggestion (+ (* 5 cg) 20000)) ;; 20000 - взято от балды
		  (left-init (- (min ul (- ug cg) (+ ug cg))
				suggestion))
		  (right-init (+ (max ul (- ug cg) (+ ug cg))
				 suggestion))
		  (root-1 (find-root left-init))
		  (root-2 (find-root right-init))
		  (X (- (sigma phi)
			(* (1- phi)
			   (Dsigma phi))))
		  (Z (* (1- phi) rg))
		  (R (+ (* Z cg cg)
			(* X phi)))
		  (c1 (* phi rl))
		  (c2 (* -2 c1 (+ ug ul)))
		  (c3 (- (* c1 (+ (expt (+ ug ul) 2)
				  (* 2 ug ul)
				  (* -1 cg cg)))
			 R))
		  (c4 (- (* 2 ug R)
			 (* 2 c1 ul
			    (- (* (+ ug ul) ug)
			       (* cg cg)))))
		  (c5 (- (* X cg cg phi)
			 (* R ug ug)
			 (* c1 ul ul (- (* cg cg)
					(* ug ug)))))
		  (coefficients (list c1 c2 c3 c4 c5))
		  (quadro-coeffs (divide-polynomial coefficients root-1 root-2)))
	     (multiple-value-bind (root-3 root-4)
		 (apply #'solve-quadratic-complex
			(mapcar (curryr #'coerce 'double-float) quadro-coeffs))
	       (when (>= *debug-level* 7)
		 (print-content t cell) (newline)
		 (format t "eigens-Afull: ~@{~d ~}~%~%" root-1 root-2 root-3 root-4))
	       (list (abs root-1) (abs root-2) (abs root-3) (abs root-4)))))))))


;; чертовски важная штука, без неё reader macro #m не заработает нормально
;; впрочем, я реализовал деление полинома сам, так как в gsll я так и не нашёл нужной функции
(setf grid:*default-grid-type* 'grid:foreign-array)

;; упрощённая функция деления полинома: производит деление столбиком.
;; coefficients - это список, начинающийся с самой большой степени.
;; деление производится на (x - divisor). Можно указать множество делителей.
;; я её написал потому что не нашёл соответствующую функцию в gcl
(defun divide-polynomial (coefficients &rest divisors)
  (cond ((emptyp divisors) coefficients)
	(t (let ((divisor (first divisors)))
	     (apply #'divide-polynomial
		    (reverse (rest (let ((x (reduce #'(lambda (acc coef)
					       ;;(format t "acc: ~a~%" acc)
					       ;;(format t "coef: ~a~%--~%" coef)
					       (if (emptyp acc)
						   (list coef)
						   (cons (+ coef (* divisor (first acc))) acc)))
					   coefficients
					   :initial-value nil)))
				     ;;(format t "!!! acc: ~a~%" x)
				     x)))
		    (rest divisors))))))

;; небольшие тесты

;;(defparameter *eigen-cell* (make-instance 'conservative-variables
;;					  :phi 0.01
;;					  :ug 25
;;					  :ul 400
;;					  :rg 1
;;					  ))
;;
;;(eigen-values-by-cell *eigen-cell* 'A)
;;(eigen-values-by-cell *eigen-cell* 'S) => -3e-7
;;(eigen-values-by-cell *eigen-cell* 'Afull)

(defun radius (cell type)
  "Спектральный радиус в методе Курганова-Тадмора"
  (let ((radius (apply #'max (mapcar #'abs (eigen-values-by-cell cell type)))))
    (when *debug-radius*
      (format t "RADIUS: ~a~%" radius))
    radius))
      

;; потоковый член в ячейке.
(defmethod flux ((vars conservative-variables))
  (let-conservative vars
    (make-instance 'conservative-variables
		   :var-vector (vector u2
				       (+ (/ (* u2 u2)
					     u1)
					  (/ (* *Cg* *Cg* u1)
					     *kappa*))
				       u4
				       (+ (/ (* u4 u4)
					     u3)
					  (/ (* *sos* *sos* u1 u3)
					     (* *kappa* (- *rl* u3)))
					  (/ (* u3 (sigma (phi vars)))
					     *rl*))))))

;; источниковый член в ячейке
(defmethod source ((vars conservative-variables))
  (let* ((d (- (ug vars) (ul vars)))
	 (sg (* -0.5 *rl* d d (signum d))))
    (make-instance 'conservative-variables
		   :var-vector (vector 0 sg 0 (- sg)))))


;; определение шага по времени

(defmethod max-timestep ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
  (if (not cu)
      *constant-timestep*
      (let ((tau 1000))
	(loop for i to (cells-max-x grid) do
	     (loop for j to (cells-max-y grid) do
		  (let* ((index (index i j))
			 (cell (cell grid index))
			 (h (max (x-size cell)
				 (y-size cell)))
			 (a+ (pseudo-vel grid index 'right))
			 (a- (pseudo-vel grid index 'left))
			 (new-tau (/ (* h cu) (max a+ a-))))
		    (when (< new-tau tau)
		      (setf tau new-tau)))))
	tau)))

(defmethod max-timestep ((grid grid2d-cv) (method (eql :kurganov-tadmor-explicit)) cu)
  (if (not cu)
      *constant-timestep*
      (let ((tau (max-timestep grid :kurganov-tadmor cu)))
	(loop for i to (cells-max-x grid) do
	     (loop for j to (cells-max-y grid) do
		  (let* ((cell (cell grid (index i j)))
			 (new-tau (/ *source-timestep-multiplier*
				     (apply #'max (eigen-values-by-cell cell 'source)))))
		    (when (< new-tau tau)
		      (setf tau new-tau)))))
	tau)))

(defun grad1-x (grid indx)
  (let* ((i (x-index indx))
	 (j (y-index indx))
	 (outborder-left (< i (min-x-index grid)))
	 (outborder-right (> i (max-x-index grid)))
	 (outborder-top (< j (min-y-index grid)))
	 (outborder-bottom (> j (max-y-index grid)))
	 (zero-gradient (make-instance 'primitive-variables
				       :number 0.0d0)))
    (cond (outborder-left zero-gradient)
	  (outborder-right zero-gradient)
	  (outborder-top zero-gradient)
	  (outborder-bottom zero-gradient)
	  (t
	   (let* ((mid-cell (cell grid indx))
		  (left-cell (unless outborder-left (cell grid (index (1- i) j))))
		  (right-cell (unless outborder-right (cell grid (index (1+ i) j))))
		  (xm (x-coord mid-cell))
		  (xl (unless outborder-left (x-coord left-cell)))
		  (xr (unless outborder-right (x-coord right-cell))))
	     (minmod 'primitive-variables
		     (pvar* *minmod-coef*
			    (pvar/ (pvar- (conservative->primitive right-cell)
					  (conservative->primitive mid-cell))
				   (- xr xm)))
		     (pvar* *minmod-coef*
			    (pvar/ (pvar- (conservative->primitive mid-cell)
					  (conservative->primitive left-cell))
				   (- xm xl)))
		     (pvar/ (pvar- (conservative->primitive right-cell)
				   (conservative->primitive mid-cell))
			    (* 2 (- xr xl)))
		     ))))))

(defmethod update-gradients ((grid grid2d-cv))
  (loop for i from (min-x-index/border grid) to (max-x-index/border grid) do
       (loop for j from (min-y-index/border grid) to (max-y-index/border grid) do
	    (let ((index (index i j)))
	      (when (valid-index-p grid index)
		(setf (gradients (cell grid index))
		      (grad1-x grid index)))))))

(defmethod copy-gradients ((oldgrid grid2d-cv) (newgrid grid2d-cv))
  (loop for i from (min-x-index/border oldgrid) to (max-x-index/border oldgrid) do
       (loop for j from (min-y-index/border oldgrid) to (max-y-index/border oldgrid) do
	    (let ((index (index i j)))
	      (when (valid-index-p oldgrid index)
		(setf (gradients (cell newgrid index))
		      (gradients (cell oldgrid index))))))))

(defun grad2-phi (grid indx)
  (let* ((i (x-index indx))
	 (j (y-index indx))
	 (mid-cell (cell grid indx))
	 (left-cell (cell grid (index (1- i) j)))
	 (right-cell (cell grid (index (1+ i) j)))
	 (xm (x-coord mid-cell))
	 (xl (x-coord left-cell))
	 (xr (x-coord right-cell))
	 (dxl (- xm xl))
	 (dxr (- xr xm))
	 (mid-phi (phi mid-cell))
	 (left-phi (phi left-cell))
	 (right-phi (phi right-cell)))
    (/ (+ (/ (* dxr (- mid-phi left-phi))
	     dxl)
	  (/ (* dxl (- right-phi mid-phi))
	     dxr))
       (+ dxr dxl))))

(defun As*dU/dx (grid index)
  (let* ((cell (cell grid index))
	 (P (P cell))
	 (grad-phi (grad2-phi grid index)))
    (cvar* P (make-instance 'conservative-variables
			    :var-vector (vector 0 grad-phi 0 (- grad-phi))))))

;; в функцию выстрела необходимо передавать градиенты, потому что
;; между предиктором и корректором мы их фиксируем
(defun shot (cell direction)
  (let ((pvars (conservative->primitive cell))
	(gx (gradients cell))
	(hx (x-size cell)))
    (let ((shot (primitive->conservative
		 (case direction
		   ((left) (pvar- pvars (pvar* 0.5 hx gx)))
		   ((right) (pvar+ pvars (pvar* 0.5 hx gx)))
		   (otherwise (error 'shot "Неизвестное направление для стрельбы"))))))
      (when *debug-shots*
	(format t "~%SHOT: ~a~%" direction)
	(format t "SHOT: from-cell:~%")
	(print-content t cell)
	(format t "~%SHOT: result:~%")
	(print-content t shot)
	(newline)
	)
      shot)))

(defun left-shot (cell)
  (shot cell 'left))

(defun right-shot (cell)
  (shot cell 'right))

;; псевдоскорость
(defun pseudo-vel (grid indx direction)
  (let* ((i (x-index indx))
	 (j (y-index indx))
	 (l-indx (index (1- i) j))
	 (r-indx (index (1+ i) j))
	 (outborder-left (< i 0))
	 (outborder-right (>= i (cells-amount-x grid)))
	 (mid-cell (cell grid indx))
	 (left-cell (unless outborder-left (cell grid l-indx)))
	 (right-cell (unless outborder-right (cell grid r-indx))))
    (when *debug-newton*
      (format t "~%NEWTON: Считаем псевдоскорость в ячейке (~d ~d)~%" i j))
    (let ((vel (case direction
		 ((right) (max (radius (primitive->conservative (right-shot mid-cell)) 'A)
			       (radius (primitive->conservative (left-shot right-cell)) 'A)
			       (radius mid-cell 'Afull)
			       (radius right-cell 'Afull)
			       ))
		 ((left) (max (radius (primitive->conservative (right-shot left-cell)) 'A)
			      (radius (primitive->conservative (left-shot mid-cell)) 'A)
			      (radius left-cell 'Afull)
			      (radius right-cell 'Afull))))))
      (when *debug-pseudo-vel* (format t "~%PSEUDO-VEL: ~d~%" vel))
      vel)))

;; потоки через границы ячейки

(defmethod border-flux ((grid grid2d-cv) (indx+ cell-index) (tau number) (border (eql :left)))
  (let* ((indx- (index (1- (x-index indx+))
		       (y-index indx+)))
	 (u+ (primitive->conservative (left-shot (cell grid indx+))))
	 (u- (primitive->conservative (right-shot (cell grid indx-))))
	 (F- (flux u+))
	 (F+ (flux u-))
	 (a (pseudo-vel grid indx- 'right)))
    (cvar* 0.5
	   (cvar- (cvar+ F+ F-)
		  (cvar* a (cvar- u+ u-))))))

(defmethod border-flux ((grid grid2d-cv) (indx cell-index) (tau number) (border (eql :right)))
  (let ((indx (index (1+ (x-index indx))
		     (y-index indx))))
    (border-flux grid indx tau :left)))

(defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number) (method (eql :kurganov-tadmor-explicit)))
  (labels
      ((predictor/corrector (grid newgrid tau)
	 (loop for j to (cells-max-y grid) do
	      (loop for i to (cells-max-x grid)
		 for dF- = (border-flux grid (index i j) tau :left) then dF+
		 for dF+ = (border-flux grid (index i j) tau :right)
		 do
		   (let* ((u (cell grid (index i j)))
			  (S (source u))
			  (hx (x-size u))
			  (dF (cvar- dF+ dF-))
			  (P*Dphi (As*dU/dx grid (index i j)))
			  (nu (cvar- (cvar+ u (cvar* tau S))
				     (cvar/ (cvar* tau dF)
					    hx)
				     (cvar* tau P*Dphi))))
		     (when *debug-predictor/corrector*
		       (format t "~%~@{~{PRED-CORR: ~a = ~d~%~}~}"
			       `((i j) (,i ,j))
			       `(oldcell-cv ,(print-conservative nil u :with-headers nil))
			       `(oldcell-pv ,(print-content nil u :with-headers nil))
			       `(tau ,tau)
			       `(dF+ ,(print-conservative nil dF+ :with-headers nil))
			       `(dF- ,(print-conservative nil dF- :with-headers nil))
			       `(dF ,(print-conservative nil dF :with-headers nil))
			       `(S ,(print-conservative nil S :with-headers nil))
			       `(P*Dphi ,(print-conservative nil P*Dphi :with-headers nil))
			       `(newcell-cv ,(print-conservative nil nu :with-headers nil))
			       `(newcell-pv ,(print-content nil nu :with-headers nil)))
		       (enbreak "Шаг предиктора/корректора"))
		     (setf (cell newgrid (index i j)) nu)
		     ))
	      (gridbg "X: predictor/corrector~%" newgrid))))
    ;; считаем градиенты исходной сетки
    (update-gradients grid)
    (gridbg "5: update-gradients~%" grid)
    ;; считаем приблизительные значения в (n+1)-й момент времени
    (dbgmsg "------------------------------------------------------------~%")
    (dbgmsg "Делаем предиктор~%")
    (when (or *debug-shots* *debug-predictor/corrector*) (print-grid grid :print-out-border nil))
    (let ((predgrid (clone grid)))
      (predictor/corrector grid predgrid tau)
      (gridbg "6: predgrid-before~%" predgrid)
      ;; в качестве сетки, полученной на шаге-предикторе, берём полусумму старой и новой сеток
      (loop for i to (cells-max-x grid) do
	   (loop for j to (cells-max-y grid) do
		(let* ((index (index i j))
		       (old-vars (primitive->conservative (cell grid index)))
		       (new-vars (primitive->conservative (cell predgrid index))))
		  (setf (cell predgrid index)
			(cvar* 0.5d0 (primitive->conservative 
				      (cvar+ old-vars new-vars)))))))
      (gridbg "7: predgrid-after~%" predgrid)
      ;; а градиенты оставляем теми же, что были на предикторе
      (copy-gradients grid predgrid)
      (gridbg "8: copy-gradients~%" predgrid)
      ;; надо удовлетворять граничным условиям
      (update-boundary-cells predgrid)
      (gridbg "8.1: update-boundary-cells~%" predgrid)
      ;; повторяем тот же шаг, но с сеткой-предиктором и старыми коэффиэнтами - это и есть корректор
      (dbgmsg "------------------------------------------------------------~%")
      (dbgmsg "Делаем корректор~%")
      (when (or *debug-shots* *debug-predictor/corrector*) (print-grid predgrid :print-out-border nil))
      (predictor/corrector predgrid newgrid tau)
      (gridbg "9: newgrid~%" newgrid)
      )
    ;;newgrid ; в будущем планирую возвращать сетку, а не писваивать значения. Может зря.
    ))

;; вычисление невязки
(defmethod count-convergence ((oldcell cell2d-cv) (newcell cell2d-cv))
  (let ((oldcell (conservative->primitive oldcell))
	(newcell (conservative->primitive newcell)))
    (apply #'max (vector->list (var-vector (pvar- newcell oldcell))))))

(defmethod count-convergence ((oldgrid grid2d-cv) (newgrid grid2d-cv))
  (setf (convergence newgrid) 0.0d0)
  (loop for j from (min-y-index oldgrid) to (max-y-index oldgrid) do
       (loop for i from (min-x-index oldgrid) to (max-x-index oldgrid)
	  for index = (index i j)
	  for oldcell = (cell oldgrid index)
	  for newcell = (cell newgrid index)
	  for conv = (count-convergence oldcell newcell)
	  do
	    (setf (convergence newcell) conv)
	    (when (> conv (convergence newgrid))
	      (setf (convergence newgrid) conv)))))

;; солвер для неявной схемы курганова-тадмора
;(defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number) (method (eql :kurganov-tadmor)))
;  (labels
;      ((predictor (grid tau)
;	 (let ((newgrid (clone grid)))
;	   (loop for j to (cells-max-y grid) do
;		(loop for i to (cells-max-x grid)
;		   for dF+ = (border-flux grid (index i j) tau :right)
;		   for dF- = (border-flux grid (index i j) tau :left) then dF+
;		   do
;		     (let* ((dF (cvar- dF+ dF-))
;			    (u (cell grid (index i j)))
;			    (hx (x-size u))
;			    (s (source u)))
;		       (labels (;; функция, для которой мы будем искать нули методом Ньютона
;				(eq-system (x1 x2 x3 x4)
;				  (let* ((x (make-instance 'conservative-variables
;							   :var-vector (vector x1 x2 x3 x4)))
;					 (result (cvar- (cvar+ (cvar/ (cvar- x u)
;								      tau)
;							       (cvar/ dF hx))
;							(cvar* 0.5 (cvar+ s (source x))))))
;				    (let-conservative result
;				      (values u1 u2 u3 u4))))
;				;; частные производеные этой функции (якобиан)
;				(eq-system-df (x1 x2 x3 x4)
;				  (let* ((D (+ (/ x2 x1)
;					       (/ x4 x3)))
;					 (F (* 0.5 *rl* D (signum D)))
;					 (F1 (/ (* F x2) (* x1 x1)))
;					 (F2 (/ F x1))
;					 (F3 (/ (* F x4) (* x3 x3)))
;					 (F4 (/ F x3))
;					 (TT (/ 1 tau)))
;				      (values TT 0.0 0.0 0.0 ;;
;					      F1 (- TT F2) (- F3) F4 ;;
;					      0.0 0.0 TT 0.0 ;;
;					      (- F1) F2 F3 (- TT F4))))
;				;; собирательная функция (необходима для GSLL)
;				(eq-system-fdf (x1 x2 x3 x4)
;				  (multiple-value-bind (v1 v2 v3 v4)
;				      (eq-system x1 x2 x3 x4)
;				    (multiple-value-bind (j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)
;					(eq-system-df x1 x2 x3 x4)
;				      (values v1 v2 v3 v4 j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)))))
;			 
;			 (let-conservative u
;			   (let ((max-iter *multi-root-iter*)
;				 (solver (make-multi-dimensional-root-solver-fdf
;					  +newton-mfdfsolver+
;					  (list #'eq-system #'eq-system-df #'eq-system-fdf)
;					  (grid:make-foreign-array
;					   'double-float :initial-contents (list u1 u2 u3 u4)))))
;			     (loop for iter from 0
;				while (if (< iter max-iter)
;					  (or (zerop iter)
;					      (not (multiroot-test-residual solver *multi-root-est-error*))))
;				do
;				  (format t "iter ~a: ~a~%" iter (vector-double-float->vector (solution solver)))
;				  (iterate solver)
;				finally
;				  ;;(format t "solution: ~a~%" (solution solver)
;				  (setf (cell newgrid (index i j))
;					      (make-instance 'conservative-variables
;							     :var-vector (vector-double-float->vector
;									  (solution solver)))))))))))
;	   newgrid))
;
;       (corrector (grid predgrid tau)
;	 (let ((midgrid (clone grid))
;	       (newgrid (clone grid)))
;	   (loop for j to (cells-max-y grid) do
;		(loop for i to (cells-max-x grid) do
;		     (setf (cell midgrid (index i j))
;			   (cvar* 0.5 (cvar+ (cell grid (index i j))
;					     (cell predgrid (index i j)))))))
;	   (loop for j to (cells-max-y grid) do
;		(loop for i to (cells-max-x grid)
;		   for index = (index i j)
;		   for dF+ = (border-flux midgrid index tau :right)
;		   for dF- = (border-flux midgrid index tau :left) then dF+
;		   do
;		     (let* ((u (cell grid index))
;			    (S (cvar* 0.5 (cvar+ (source u)
;						 (source (cell predgrid index)))))
;			    (hx (x-size u))
;			    (dF (cvar- dF+ dF-)))
;		       (setf (cell newgrid index)
;			     (cvar- (cvar+ u (cvar* tau S))
;				    (cvar/ (cvar* tau dF)
;					   hx))))))
;	   newgrid)))
;    (let ((predgrid (predictor grid tau)))
;      (setf newgrid predgrid)
;      ;;(corrector grid predgrid tau)
;
;      (format t "PREDICTOR~%~%")
;      (print-grid newgrid :print-out-border t)
;      (newline)(newline)
;      )))
       

				      
			      

;(defmacro list-of-vars (prefix argn &rest body)
;  (let* ((fmt (format nil "~a~~d" (string-upcase (any->string prefix))))
;	 (symbols (mapcar #'intern (loop for i from 1 to (eval argn)
;				      collect (format nil fmt i)))))
;    symbols))
;
;(defun inc-vector-elem (vec i dx)
;  (let ((vec (copy-seq vec)))
;    (incf (aref vec i) dx)
;    vec))
;
;(defun make-jacobian (func &rest argv)
;  (let ((tt 1d-5)
;	(argc (1- (length argv))))
;    (format t "1~%")
;    (multiple-value-bind ((list-of-vars "f" argc) (apply func argv))
;	(loop for i to argc do
;	     (append (loop for j to argc do
;			  (collect (/ (- (apply func (inc-vector-elem argv j tt)) (apply func argv))
;				      tt))))))))
;