Browse code

Тест на монотонном поле + много отладки

Слот градиентов в класс консервативных переменных.
Валидация индексов: предикаты местоположения индекса относительно сетки.
Стандартизированные функции для определения минимальных и максимальных индексов сеток.
Включил в зависимости систему cl-ana и memoization - залог на будущее. Пока их можно свободно выбросить.
Вывод сообщений для дебага.
Функции копирования и обновления градиентов.
Фиксация градиентов на корректоре в методе Курганова-Тадмора.
Метод Курганова-Тадмора теперь явный.

Dmitrii Kashin authored on 18/05/2015 14:49:20
Showing 13 changed files
... ...
@@ -62,7 +62,7 @@
62 62
 ;; -------------------- ДЛЯ КОНСЕРВАТИВНЫХ ПЕРЕМЕННЫХ
63 63
 
64 64
 (defmethod update-boundary-cells ((grid grid2d-cv))
65
-  (loop for i from 0 to (cells-max-x grid) do
65
+  (loop for i from (min-x-index grid) to (max-x-index grid) do
66 66
        (funcall (symbol-function (bottom-boundary-condition grid))
67 67
 		grid
68 68
 		(index i -1)
... ...
@@ -85,11 +85,13 @@
85 85
   (let ((i (x-index index))
86 86
 	(j (y-index index)))
87 87
     (setf (cell grid index)
88
-	  (cell grid (case border
89
-		       ((left) (index (1+ i) j))
90
-		       ((right) (index (1- i) j))
91
-		       ((bottom) (index i (1+ j)))
92
-		       ((top) (index i (1- j))))))))
88
+	  (let ((cell (cell grid (case border
89
+				   ((left) (index (1+ i) j))
90
+				   ((right) (index (1- i) j))
91
+				   ((bottom) (index i (1+ j)))
92
+				   ((top) (index i (1- j)))))))
93
+	    (make-instance 'conservative-variables
94
+			   :var-vector (var-vector cell))))))
93 95
 
94 96
 (defmethod fixed ((grid grid2d-cv) (index cell-index) border)
95 97
   (noop))
... ...
@@ -26,12 +26,44 @@
26 26
   ((dimensions :initarg :dimensions
27 27
 	       :initform 2
28 28
 	       :accessor dimensions)
29
-   (headers/cell2d-cv :initform '()
29
+   ;; gradients содержит массив примитивных переменных и их градиентов
30
+   (gradients :initform (make-dummy-variables 'primitive-variables)
31
+	      :initarg :gradients
32
+	      :accessor gradients)
33
+   (headers/cell2d-cv :initform '(grad-P grad-phi grad-ug grad-ul)
30 34
 		      :allocation :class)))
31 35
 
36
+;; акцессоры и мутаторы градиентов (а нужны ли вообще? ну, пускай будут)
37
+(defun grad-phi (cell)
38
+  (phi (gradients cell)))
39
+
40
+(defun grad-P (cell)
41
+  (P (gradients cell)))
42
+
43
+(defun grad-ug (cell)
44
+  (ug (gradients cell)))
45
+
46
+(defun grad-ul (cell)
47
+  (ul (gradients cell)))
48
+
49
+(defun (setf grad-phi) (new cell)
50
+  (setf (phi (gradients cell)) new))
51
+
52
+(defun (setf grad-P) (new cell)
53
+  (setf (P (gradients cell)) new))
54
+
55
+(defun (setf grad-ug) (new cell)
56
+  (setf (ug (gradients cell)) new))
57
+
58
+(defun (setf grad-ul) (new cell)
59
+  (setf (ul (gradients cell)) new))
60
+
61
+(make-header-for-class cell2d-cv)
62
+
32 63
 (print-functions-for-class cell2d-cv 
33 64
 			   ((inherit . cell-coordinates)
34 65
 			    (inherit . conservative-variables)
66
+			    (inherit . cell2d-cv)
35 67
 			    (inherit . cell-size)
36 68
 			    (inherit . cell-base)))
37 69
 
... ...
@@ -39,6 +71,8 @@
39 39
   (make-instance type
40 40
 		 :coord-vector (vector 'undef 'undef)
41 41
 		 :size-vector (vector 'undef 'undef)
42
+		 :gradients (make-instance 'primitive-variables
43
+					   :var-vector (vector 'undef 'undef 'undef 'undef))
42 44
 		 :number 'undef))
43 45
 
44 46
 (defmethod clone ((cell cell2d-cv))
... ...
@@ -50,3 +84,11 @@
50 50
 		 :ug (ug cell)
51 51
 		 :ul (ul cell)
52 52
 		 ))
53
+
54
+
55
+;; строго говоря, не правильно добавлять градиенты в cell2d-cv. Надо
56
+;; было определить класс cell1d-cv-order2, добавить в него градиенты
57
+;; по x и наследовать от cell1d-cv, а затем cell2d-cv-order2, и
58
+;; наследовать его от класса cell1d-cv-order2 (а также от cell2d-cv)
59
+
60
+;; но я сейчас быстро хакаю - магистерская на носу
... ...
@@ -22,7 +22,11 @@
22 22
 (defparameter *newton-max-iter* 50)
23 23
 (defparameter *newton-est-error* 1.0d-8)
24 24
 
25
-(defparameter *multi-root-iter* 1000)
25
+(defparameter *multi-root-iter* 10)
26 26
 (defparameter *multi-root-est-error* 1.0d-8)
27 27
 
28 28
 (defparameter *constant-timestep* .5d-8)
29
+
30
+(defparameter *source-timestep-multiplier* 0.2d0)
31
+
32
+(defparameter *debug-level* 5)
... ...
@@ -283,7 +283,7 @@
283 283
 метода `method'. Метод может быть одним из следующих
284 284
 значений: :lax-frederichs :kurganov-tadmor"))
285 285
 
286
-(defgeneric make-step (grid newgrid tau)
286
+(defgeneric make-step (grid newgrid tau method)
287 287
   (:documentation "Делает в сетке `grid' шаг по времени величиной
288 288
   `tau'. Результат записывается в `newgrid'"))
289 289
 
... ...
@@ -40,42 +40,120 @@
40 40
   ;; создать сам массив
41 41
   (setf (cells grid)
42 42
 	(make-array (list (2+ (cells-amount-x grid))
43
-			  (2+ (cells-amount-y grid)))
44
-		    :element-type 'cell2d-cv))
45
-  ;; создать все его элементы
46
-  (loop for i from -1 to (1+ (cells-max-x grid)) do
47
-       (loop for j from -1 to (1+ (cells-max-y grid)) do
48
-	    ;;(format t "~d:~d" i j) (newline)
49
-	    (setf (cell grid (index i j))
50
-		  (make-dummy-cell 'cell2d-cv))
51
-	    ))
52
-  )
43
+			  (2+ (cells-amount-y grid)))))
44
+		    ;;:element-type 'cell2d-cv))
45
+  ;; создать все его элементы и их градиенты
46
+  (let ((xi-min -1)
47
+	(yi-min -1)
48
+	(xi-max (cells-amount-x grid))
49
+	(yi-max (cells-amount-y grid)))
50
+    (loop for i from (min-x-index/border grid) to (max-x-index/border grid) do
51
+	 (loop for j from (min-y-index/border grid) to (max-y-index/border grid) do
52
+	      (let ((index (index i j)))
53
+		(setf (cell grid index)
54
+		      (if (valid-index-p grid index)
55
+			  (make-dummy-cell 'cell2d-cv)
56
+			  'outgrid)))))))
57
+
58
+;; проверка, находится ли индекс в допустимых диапазонах
59
+(defun grid-index-type (grid index)
60
+  (let ((xi (x-index index))
61
+	(yi (y-index index))
62
+	(xi-min -1)
63
+	(yi-min -1)
64
+	(xi-max (cells-amount-x grid))
65
+	(yi-max (cells-amount-y grid))
66
+	)
67
+    (cond ((and (<= 0 xi (1- xi-max))
68
+		(<= 0 yi (1- yi-max)))
69
+	   'internal)
70
+	  ((or (and (or (equal xi xi-min)
71
+			(equal xi xi-max))
72
+		    (<= 0 yi (1- yi-max)))
73
+	       (and (or (equal yi yi-min)
74
+			(equal yi yi-max))
75
+		    (<= 0 xi (1- xi-max))))
76
+	   'border)
77
+	  (t 'outgrid))))
78
+
79
+(defun border-index-p (grid index)
80
+  (equal (grid-index-type grid index) 'border))
81
+
82
+(defun internal-index-p (grid index)
83
+  (equal (grid-index-type grid index) 'internal))
84
+
85
+(defun outgrid-index-p (grid index)
86
+  (equal (grid-index-type grid index) 'outgrid))
87
+
88
+(defun valid-index-p (grid index)
89
+  (not (outgrid-index-p grid index)))
90
+
91
+;; выдаёт ошибку, если индекс вне допустимых диапазонов
92
+(defun index-validation (grid index)
93
+  (unless (valid-index-p grid index)
94
+    (error "Индекс за пределами сетки: (index ~a ~a)~%" (x-index index) (y-index index))))
53 95
 
54 96
 ;; Акцессоры и мутаторы для элементов сетки
55 97
 ;; note: у меня такое чувство, будто мутаторами я пользоваться не буду
56 98
 
57 99
 (defmethod cell ((grid grid2d-cv) (index cell-index))
58
-  (unless (and (<= -1 (x-index index) (cells-amount-x grid))
59
-	       (<= -1 (y-index index) (cells-amount-y grid)))
60
-    (error "Индекс за пределами сетки: (index ~a ~a)~%" (x-index index) (y-index index)))
100
+  (index-validation grid index)
61 101
   (aref (cells grid) (1+ (x-index index)) (1+ (y-index index))))
62 102
 
63 103
 (defmethod (setf cell) ((new cell2d-cv) (grid grid2d-cv) (index cell-index))
104
+  (index-validation grid index)
64 105
   (setf (aref (cells grid) (1+ (x-index index)) (1+ (y-index index)))
65 106
 	(clone new)))
66 107
 
67 108
 (defmethod (setf cell) ((new conservative-variables) (grid grid2d-cv) (index cell-index))
109
+  (index-validation grid index)
68 110
   (setf (var-vector (cell grid index))
69 111
 	(var-vector new)))
70 112
 
71 113
 (defmethod (setf cell) ((new cell-coordinates) (grid grid2d-cv) (index cell-index))
114
+  (index-validation grid index)
72 115
   (setf (coord-vector (cell grid index))
73 116
 	(coord-vector new)))
74 117
 
75 118
 (defmethod (setf cell) ((new cell-size) (grid grid2d-cv) (index cell-index))
119
+  (index-validation grid index)
76 120
   (setf (size-vector (cell grid index))
77 121
 	(size-vector new)))
78 122
 
123
+;; Этот метод используется для задания символов 'OUTGRID в угловых
124
+;; ячейках сетки, поэтому валидация не нужна.
125
+(defmethod (setf cell) ((new symbol) (grid grid2d-cv) (index cell-index))
126
+  (setf (aref (cells grid) (1+ (x-index index)) (1+ (y-index index))) new))
127
+	
128
+  
129
+ 
130
+
131
+;; дополнительные акцессоры для сетки
132
+
133
+(defmethod min-x-index ((grid grid2d-cv))
134
+  0)
135
+
136
+(defmethod max-x-index ((grid grid2d-cv))
137
+  (1- (cells-amount-x grid)))
138
+
139
+(defmethod min-x-index/border ((grid grid2d-cv))
140
+  (1- (min-x-index grid)))
141
+
142
+(defmethod max-x-index/border ((grid grid2d-cv))
143
+  (1+ (max-x-index grid)))
144
+
145
+(defmethod min-y-index ((grid grid2d-cv))
146
+  0)
147
+
148
+(defmethod max-y-index ((grid grid2d-cv))
149
+  (1- (cells-amount-y grid)))
150
+
151
+(defmethod min-y-index/border ((grid grid2d-cv))
152
+  (1- (min-y-index grid)))
153
+
154
+(defmethod max-y-index/border ((grid grid2d-cv))
155
+  (1+ (max-y-index grid)))
156
+
79 157
 ;; дополнительные акцессоры используются из grid2d-primitive.lisp
80 158
 ;; функции инициализации я подцепляю оттуда же
81 159
 
... ...
@@ -203,24 +281,15 @@
203 203
   (newline)
204 204
   ;; ячейки
205 205
   (loop for i from -1 to (1+ (cells-max-x grid)) do
206
-;  (loop for i from -1 to (1+ max-x) do ; bit-hack для тестов
207 206
        (loop for j from -1 to (1+ (cells-max-y grid)) do
208
-	    (let ((cell (cell grid (index i j)))
207
+	    (let ((index (index i j))
209 208
 		  (print-flag nil))
210
-	      ;; исключаем углы
211
-	      (unless (and (not (<= 0 i (cells-max-x grid)))
212
-			   (not (<= 0 j (cells-max-y grid))))
213
-		(cond (;; границы
214
-		       (or (= i -1)
215
-			   (= i (1+ (cells-max-x grid)))
216
-			   (= j -1)
217
-			   (= j (1+ (cells-max-y grid))))
218
-		       (when print-out-border (setf print-flag t)))
219
-		      (t ;; центральные
220
-		       (setf print-flag t))))
209
+	      (case (grid-index-type grid index)
210
+		((border) (when print-out-border (setf print-flag t)))
211
+		((internal) (setf print-flag t)))
221 212
 	      (when print-flag
222 213
 		(format t "~d ~d " i j)
223
-	        (print-content t cell :with-headers nil)
214
+	        (print-content t (cell grid (index i j)) :with-headers nil)
224 215
 		(newline))))))
225 216
 
226 217
 (defmethod clone ((grid grid2d-cv))
... ...
@@ -88,6 +88,8 @@
88 88
 (defun nodes-max-y (grid)
89 89
   (cells-amount-y grid))
90 90
 
91
+
92
+
91 93
 ;; функции, необходимые для инициализации сетки
92 94
 ;; сразу скажу, что это старые функции, которые выполняются лишь один раз, и больше не требуются
93 95
 ;; вообще, если по-хорошему, то надо заменить все хэши на двумерные массивы, но мне лень
... ...
@@ -17,10 +17,13 @@
17 17
 ;(push "/home/freehck/prog/maxima-5.35.1/src/" asdf:*central-registry*)
18 18
 ;(asdf:load-system :maxima)
19 19
 
20
+(ql:quickload "cl-ana")
21
+
22
+
20 23
 (asdf::load-asd "clf.asd")
21 24
 ;(asdf::load-asd "alexandria.asd")
22 25
 
23 26
 (asdf:load-system :alexandria)
24 27
 (asdf:load-system :gsll)
28
+(asdf:load-system :cl-ana)
25 29
 (asdf:load-system :clf)
26
-
... ...
@@ -2,8 +2,7 @@
2 2
 
3 3
 (defun march (grid &key output-every-step output-every-time finish-time finish-step
4 4
 		     print-out-border suppress-export
5
-		     (max-timestep #'max-timestep)
6
-		     (make-step #'make-step)
5
+		     (method :kurganov-tadmor-explicit)
7 6
 		     (cu nil))
8 7
   ;; Проверка начальных условий и вывод информации о процессе
9 8
   (format t "Маршевая функция запущена.~%")
... ...
@@ -26,56 +25,63 @@
26 26
     ;; если finish-step установлен в 0, то больше ничего не делаем
27 27
     (when (<= finish-step 0)
28 28
       (setf stop-flag t))
29
-    (loop until stop-flag
30
-       do (let* ((print-flag nil)
31
-		 (time (grid-time grid))
32
-		 (tau (funcall max-timestep grid :kurganov-tadmor cu))
33
-		 (iterations (iteration grid)))
34
-	    ;; изменяем шаг по времени, если надо делать вывод в
35
-	    ;; определённый момент
36
-	    (when output-every-time
37
-	      ;; гарантируем, что не перескочим сразу несколько
38
-	      ;; моментов для вывода по времени
39
-	      (when (> tau output-every-time)
40
-		(setf tau output-every-time))
41
-	      ;; условие перехода через момент времени, когда мы
42
-	      ;; обязаны сделать вывод
43
-	      (when (< (mod (+ time tau) output-every-time)
44
-		       (mod time output-every-time))
45
-		(decf tau (mod (+ time tau) output-every-time))
46
-		(setf print-flag t))
47
-	      )
48
-	    ;; изменяем шаг по времени, если близки ко времени
49
-	    ;; окончанию расчёта
50
-	    (when (and finish-time
51
-		       (> (+ time tau) finish-time))
52
-	      (setf tau (- finish-time time))
53
-	      (setf print-flag t)
54
-	      (setf stop-flag t))
55
-	    ;; проверяем, не надо ли делать вывод после этой итерации
56
-	    (when (and output-every-step
57
-		       (zerop (mod (1+ iterations) output-every-step)))
58
-	      (setf print-flag t))
59
-	    ;; проверяем, является ли следующий шаг последним
60
-	    (when (>= (1+ iterations) finish-step)
61
-	      (setf print-flag t)
62
-	      (setf stop-flag t))
29
+    (loop until stop-flag do
30
+	 (update-gradients grid)
31
+	 (let* ((print-flag nil)
32
+		(time (grid-time grid))
33
+		(tau (max-timestep grid method cu))
34
+		(iterations (iteration grid)))
35
+	   ;; изменяем шаг по времени, если надо делать вывод в
36
+	   ;; определённый момент
37
+	   (when output-every-time
38
+	     ;; гарантируем, что не перескочим сразу несколько
39
+	     ;; моментов для вывода по времени
40
+	     (when (> tau output-every-time)
41
+	       (setf tau output-every-time))
42
+	     ;; условие перехода через момент времени, когда мы
43
+	     ;; обязаны сделать вывод
44
+	     (when (< (mod (+ time tau) output-every-time)
45
+		      (mod time output-every-time))
46
+	       (decf tau (mod (+ time tau) output-every-time))
47
+	       (setf print-flag t))
48
+	     )
49
+	   ;; изменяем шаг по времени, если близки ко времени
50
+	   ;; окончанию расчёта
51
+	   (when (and finish-time
52
+		      (> (+ time tau) finish-time))
53
+	     (setf tau (- finish-time time))
54
+	     (setf print-flag t)
55
+	     (setf stop-flag t))
56
+	   ;; проверяем, не надо ли делать вывод после этой итерации
57
+	   (when (and output-every-step
58
+		      (zerop (mod (1+ iterations) output-every-step)))
59
+	     (setf print-flag t))
60
+	   ;; проверяем, является ли следующий шаг последним
61
+	   (when (>= (1+ iterations) finish-step)
62
+	     (setf print-flag t)
63
+	     (setf stop-flag t))
64
+	   
65
+	   ;; шаг
66
+	   (format t "Iter: ~d, Tau: ~d, Time: ~d~%" iterations tau time)
67
+	   (make-step grid newgrid tau method)
68
+	   ;; заботимся о выставлении параметров сетки
69
+	   (incf (iteration newgrid))
70
+	   (incf (grid-time newgrid) tau)
71
+	   ;; печать
72
+	   (when (and (not suppress-export)
73
+		      print-flag)
74
+	     (export-grid newgrid :print-out-border print-out-border))
75
+	   
76
+	   ;; меняем сетки местами
77
+	   (setf grid (clone newgrid))
63 78
 
64
-	    ;; шаг
65
-	    (funcall make-step grid newgrid tau)
66
-	    ;; заботимся о выставлении параметров сетки
67
-	    (incf (iteration newgrid))
68
-	    (incf (grid-time newgrid) tau)
69
-	    ;; печать
70
-	    (when (and (not suppress-export)
71
-		       print-flag)
72
-	      (export-grid newgrid :print-out-border print-out-border))
73
-
74
-	    ;; меняем сетки местами
75
-	    (setf grid (clone newgrid))
76
-	    ;(update-gradients-conv grid)
77
-	    ))
79
+	   ))
78 80
 
81
+    (when (>= *debug-level* 10)
82
+	     (format t "10: final grid~%")
83
+	     (print-grid newgrid :print-out-border t)
84
+	     (format t "----~%~%"))
85
+    
79 86
     ;; по завершении маршевой функции возвращаем последнюю полученную сетку
80 87
     grid
81 88
     ))
... ...
@@ -3,10 +3,13 @@
3 3
 (defpackage :clf
4 4
   (:use :common-lisp
5 5
 	:alexandria
6
-	:gsll)
6
+	:gsll
7
+	;;:memoization ; в будущем я смогу использовать эту подсистему cl-ana, чтобы ускорить работу методов
8
+	)
7 9
   (:shadow
8 10
    ;; следующие символы конфликтовали между alexandria и gsll
9 11
    mean variance median standard-deviation factorial
10 12
    ;; следующие символы конфликтовали с теми, что я определяю в этом пакете
11 13
    sigma))
12 14
 
15
+
... ...
@@ -24,17 +24,28 @@
24 24
 			       (* phi (- (sigma phi)
25 25
 					 rl ul ul)))
26 26
 			    (* phi rl)))))
27
-	   (list (abs (- ug cg))
28
-		 (abs (+ ug cg))
29
-		 (abs (+ ul cl))
30
-		 (abs (- ul cl)))))
27
+	   (let ((x1 (- ug cg))
28
+		 (x2 (+ ug cg))
29
+		 (x3 (+ ul cl))
30
+		 (x4 (- ul cl)))
31
+	     (when (>= *debug-level* 7)
32
+	       (print-content t cell) (newline)
33
+	       (format t "eigens-A: ~@{~d ~}~%~%" x1 x2 x3 x4))
34
+	     (list (abs x1) (abs x2) (abs x3) (abs x4)))))
31 35
 
32 36
 	((S source) ;; собственные числа (одно) якобиана источниковых членов
33
-	 (list (- (* (abs (- ug ul))
34
-		     (+ (/ rl
35
-			   (* phi rg))
36
-			(/ 1
37
-			   (- 1 phi)))))))
37
+	 (let* ((ug-ul (if (= ug ul)
38
+			   1.d-8
39
+			   (- ug ul)))
40
+		(x1 (- (* (abs ug-ul)
41
+			  (+ (/ rl
42
+				(* phi rg))
43
+			     (/ 1
44
+				(- 1 phi)))))))
45
+	   (when (>= *debug-level* 6)
46
+	     (print-content t cell) (newline)
47
+	     (format t "eigens-S: ~d~%~%" x1))
48
+	   (list (abs x1))))
38 49
 
39 50
 	((Afull) ;; собственные числа якобиана потоковых членов, объединённых с
40 51
 		 ;; источниками, обусловленными сужением трубки тока
... ...
@@ -109,7 +120,10 @@
109 109
 	     (multiple-value-bind (root-3 root-4)
110 110
 		 (apply #'solve-quadratic-complex
111 111
 			(mapcar (curryr #'coerce 'double-float) quadro-coeffs))
112
-	       (list root-1 root-2 root-3 root-4))))))))
112
+	       (when (>= *debug-level* 7)
113
+		 (print-content t cell) (newline)
114
+		 (format t "eigens-Afull: ~@{~d ~}~%~%" root-1 root-2 root-3 root-4))
115
+	       (list (abs root-1) (abs root-2) (abs root-3) (abs root-4)))))))))
113 116
 
114 117
 
115 118
 ;; чертовски важная штука, без неё reader macro #m не заработает нормально
... ...
@@ -149,6 +163,10 @@
149 149
 ;;(eigen-values-by-cell *eigen-cell* 'S) => -3e-7
150 150
 ;;(eigen-values-by-cell *eigen-cell* 'Afull)
151 151
 
152
+;; спектральный радиус
153
+(defun radius (cell type)
154
+  (apply #'max (mapcar #'abs (eigen-values-by-cell cell type))))
155
+
152 156
 ;; потоковый член в ячейке.
153 157
 (defmethod flux ((vars conservative-variables))
154 158
   (let-conservative vars
... ...
@@ -173,93 +191,126 @@
173 173
     (make-instance 'conservative-variables
174 174
 		   :var-vector (vector 0 sg 0 (- sg)))))
175 175
 
176
-;; градиенты
177
-(defun grad-x (grid indx order)
176
+
177
+;; определение шага по времени
178
+
179
+(defmethod max-timestep ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
180
+  (if (not cu)
181
+      *constant-timestep*
182
+      (let ((tau 1000))
183
+	(loop for i to (cells-max-x grid) do
184
+	     (loop for j to (cells-max-y grid) do
185
+		  (let* ((index (index i j))
186
+			 (cell (cell grid index))
187
+			 (h (max (x-size cell)
188
+				 (y-size cell)))
189
+			 (a+ (pseudo-vel grid index 'right))
190
+			 (a- (pseudo-vel grid index 'left))
191
+			 (new-tau (/ (* h cu) (max a+ a-))))
192
+		    (when (< new-tau tau)
193
+		      (setf tau new-tau)))))
194
+	tau)))
195
+
196
+(defmethod max-timestep ((grid grid2d-cv) (method (eql :kurganov-tadmor-explicit)) cu)
197
+  (if (not cu)
198
+      *constant-timestep*
199
+      (let ((tau (max-timestep grid :kurganov-tadmor cu)))
200
+	(loop for i to (cells-max-x grid) do
201
+	     (loop for j to (cells-max-y grid) do
202
+		  (let* ((cell (cell grid (index i j)))
203
+			 (new-tau (/ *source-timestep-multiplier*
204
+				     (apply #'max (eigen-values-by-cell cell 'source)))))
205
+		    (when (< new-tau tau)
206
+		      (setf tau new-tau)))))
207
+	tau)))
208
+
209
+(defun grad1-x (grid indx)
178 210
   (let* ((i (x-index indx))
179 211
 	 (j (y-index indx))
180
-	 (outborder-left (< i 0))
181
-	 (outborder-right (>= i (cells-amount-x grid)))
182
-	 (mid-cell (cell grid indx))
183
-	 (left-cell (unless outborder-left
184
-		      (cell grid (index (1- i) j))))
185
-	 (right-cell (unless outborder-right
186
-		       (cell grid (index (1+ i) j))))
187
-	 (DR (unless outborder-right
188
-	       (- (x-coord right-cell) (x-coord mid-cell))))
189
-	 (DL (unless outborder-left
190
-	       (- (x-coord mid-cell) (x-coord left-cell))))
191
-	 (DUR (unless outborder-right
192
-		(cvar- right-cell mid-cell)))
193
-	 (DUL (unless outborder-left
194
-		(cvar- mid-cell left-cell))))
195
-    (case order
196
-      ;; 1й порядок точности может потребоваться в заграничных ячейках
197
-      ((1) (cond (outborder-left (cvar/ DUR DR))
198
-		 (outborder-right (cvar/ DUL DL))
199
-		 (t (minmod (cvar/ DUL DL)
200
-			    (cvar/ DUR DR)))))
201
-      ;; 2й порядок может потребоваться только в центральный, поэтому cond не нужен
202
-      ((2) (cvar/ (cvar+ (cvar/ (cvar* DUL DR)
203
-				DL)
204
-			 (cvar/ (cvar* DUR DL)
205
-				DR))
206
-		  (cvar+ DL DR)))
207
-      (otherwise (error "`grad' умеет считать производные лишь до 2го порядка включительно")))))
212
+	 (outborder-left (< i (min-x-index grid)))
213
+	 (outborder-right (> i (max-x-index grid)))
214
+	 (outborder-top (< j (min-y-index grid)))
215
+	 (outborder-bottom (> j (max-y-index grid)))
216
+	 (zero-gradient (make-instance 'primitive-variables
217
+				       :number 0.0d0)))
218
+    (cond (outborder-left zero-gradient)
219
+	  (outborder-right zero-gradient)
220
+	  (outborder-top zero-gradient)
221
+	  (outborder-bottom zero-gradient)
222
+	  (t
223
+	   (let* ((mid-cell (cell grid indx))
224
+		  (left-cell (unless outborder-left (cell grid (index (1- i) j))))
225
+		  (right-cell (unless outborder-right (cell grid (index (1+ i) j))))
226
+		  (xm (x-coord mid-cell))
227
+		  (xl (unless outborder-left (x-coord left-cell)))
228
+		  (xr (unless outborder-right (x-coord right-cell))))
229
+	     (minmod (pvar/ (pvar- (conservative->primitive right-cell)
230
+				   (conservative->primitive mid-cell))
231
+			    (- xr xm))
232
+		     (pvar/ (pvar- (conservative->primitive mid-cell)
233
+				   (conservative->primitive left-cell))
234
+			    (- xm xl))))))))
208 235
 
209
-(defun grad1-x (grid indx)
210
-  (grad-x grid indx 1))
236
+(defmethod update-gradients ((grid grid2d-cv))
237
+  (loop for i from (min-x-index/border grid) to (max-x-index/border grid) do
238
+       (loop for j from (min-y-index/border grid) to (max-y-index/border grid) do
239
+	    (let ((index (index i j)))
240
+	      (when (valid-index-p grid index)
241
+		(setf (gradients (cell grid index))
242
+		      (grad1-x grid index)))))))
211 243
 
212
-(defun grad2-x (grid indx)
213
-  (grad-x grid indx 2))
244
+(defmethod copy-gradients ((oldgrid grid2d-cv) (newgrid grid2d-cv))
245
+  (loop for i from (min-x-index/border oldgrid) to (max-x-index/border oldgrid) do
246
+       (loop for j from (min-y-index/border oldgrid) to (max-y-index/border oldgrid) do
247
+	    (let ((index (index i j)))
248
+	      (when (valid-index-p oldgrid index)
249
+		(setf (gradients (cell newgrid index))
250
+		      (gradients (cell oldgrid index))))))))
214 251
 
215
-;; источниковый член сужения трубки тока за счёт межфазного обмена
216
-(defun ASF (grid indx)
217
-  (let* ((cell (cell grid indx))
218
-	 (phi (phi cell))
219
-	 (rg (rg cell))
220
-	 (ug (ug cell))
221
-	 (ul (ul cell))
222
-	 (C (/ *Cg* (sqrt *kappa*)))
223
-	 (K (/ (* C C rg)
224
-	       *rl*))
225
-	 (N (/ (1- phi)
226
-	       phi))
227
-	 (grads (grad2-x grid indx)))
228
-    (let-conservative grads
229
-      (make-instance 'conservative-variables
230
-		     :var-vector (vector u3
231
-					 (+ (* (- (* C C)
232
-						  (* ug ug))
233
-					       u1)
234
-					    (* 2 ug u2)
235
-					    (* K u3))
236
-					 u4
237
-					 (+ (* C C N u1)
238
-					    (* (- (* K N)
239
-						  (* ul ul))
240
-					       u3)
241
-					    (* 2 ul u4)))))))
242
-
243
-;; todo: а правильно ли я стреляю? линейны-то примитивные переменные
244
-;; выстрелы на границу
245
-(defun shot (grid indx direction)
246
-  (let* ((cell (cell grid indx))
247
-	 (hx (x-size cell))
248
-	 (gx (grad1-x grid indx)))
249
-    (case direction
250
-      ((left) (cvar- cell (cvar* 0.5 hx gx)))
251
-      ((right) (cvar+ cell (cvar* 0.5 hx gx)))
252
-      (otherwise (error 'shot "Неизвестное направление для стрельбы")))))
252
+(defun grad2-phi (grid indx)
253
+  (let* ((i (x-index indx))
254
+	 (j (y-index indx))
255
+	 (mid-cell (cell grid indx))
256
+	 (left-cell (cell grid (index (1- i) j)))
257
+	 (right-cell (cell grid (index (1+ i) j)))
258
+	 (xm (x-coord mid-cell))
259
+	 (xl (x-coord left-cell))
260
+	 (xr (x-coord right-cell))
261
+	 (dxl (- xm xl))
262
+	 (dxr (- xr xm))
263
+	 (mid-phi (phi mid-cell))
264
+	 (left-phi (phi left-cell))
265
+	 (right-phi (phi right-cell)))
266
+    (/ (+ (/ (* dxr (- mid-phi left-phi))
267
+	     dxl)
268
+	  (/ (* dxl (- right-phi mid-phi))
269
+	     dxr))
270
+       (+ dxr dxl))))
253 271
 
254
-(defun left-shot (grid indx)
255
-  (shot grid indx 'left))
272
+(defun As*dU/dx (grid index)
273
+  (let* ((cell (cell grid index))
274
+	 (P (P cell))
275
+	 (grad-phi (grad2-phi grid index)))
276
+    (cvar* P (make-instance 'conservative-variables
277
+			    :var-vector (vector 0 grad-phi 0 (- grad-phi))))))
256 278
 
257
-(defun right-shot (grid indx)
258
-  (shot grid indx 'right))
279
+;; в функцию выстрела необходимо передавать градиенты, потому что
280
+;; между предиктором и корректором мы их фиксируем
281
+(defun shot (cell direction)
282
+  (let ((pvars (conservative->primitive cell))
283
+	(gx (gradients cell))
284
+	(hx (x-size cell)))
285
+    (primitive->conservative
286
+     (case direction
287
+       ((left) (pvar- pvars (pvar* 0.5 hx gx)))
288
+       ((right) (pvar+ pvars (pvar* 0.5 hx gx)))
289
+       (otherwise (error 'shot "Неизвестное направление для стрельбы"))))))
259 290
 
260
-;; спектральный радиус
261
-(defun radius (cell type)
262
-  (apply #'max (mapcar #'abs (eigen-values-by-cell cell type))))
291
+(defun left-shot (cell)
292
+  (shot cell 'left))
293
+
294
+(defun right-shot (cell)
295
+  (shot cell 'right))
263 296
 
264 297
 ;; псевдоскорость
265 298
 (defun pseudo-vel (grid indx direction)
... ...
@@ -270,47 +321,25 @@
270 270
 	 (outborder-left (< i 0))
271 271
 	 (outborder-right (>= i (cells-amount-x grid)))
272 272
 	 (mid-cell (cell grid indx))
273
-	 (left-cell (unless outborder-left
274
-		      (cell grid l-indx)))
275
-	 (right-cell (unless outborder-right
276
-		       (cell grid r-indx))))
273
+	 (left-cell (unless outborder-left (cell grid l-indx)))
274
+	 (right-cell (unless outborder-right (cell grid r-indx))))
277 275
     (case direction
278
-      ((right) (max (radius (right-shot grid indx) 'A)
279
-		    (radius (left-shot grid r-indx) 'A)
276
+      ((right) (max (radius (primitive->conservative (right-shot mid-cell)) 'A)
277
+		    (radius (primitive->conservative (left-shot right-cell)) 'A)
280 278
 		    (radius mid-cell 'Afull)
281 279
 		    (radius right-cell 'Afull)))
282
-      ((left) (max (radius (right-shot grid l-indx) 'A)
283
-		   (radius (left-shot grid indx) 'A)
280
+      ((left) (max (radius (primitive->conservative (right-shot left-cell)) 'A)
281
+		   (radius (primitive->conservative (left-shot mid-cell)) 'A)
284 282
 		   (radius left-cell 'Afull)
285 283
 		   (radius right-cell 'Afull))))))
286 284
 
287
-
288
-;; определение шага по времени
289
-
290
-(defmethod max-timestep ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
291
-  (if (not cu)
292
-      *constant-timestep*
293
-      (let ((tau 1000))
294
-	(loop for i to (cells-max-x grid) do
295
-	     (loop for j to (cells-max-y grid) do
296
-		  (let* ((index (index i j))
297
-			 (cell (cell grid index))
298
-			 (h (max (x-size cell)
299
-				 (y-size cell)))
300
-			 (a+ (pseudo-vel grid index 'right))
301
-			 (a- (pseudo-vel grid index 'left))
302
-			 (new-tau (/ (* h cu) (max a+ a-))))
303
-		    (when (< new-tau tau)
304
-		      (setf tau new-tau)))))
305
-	tau)))
306
-
307 285
 ;; потоки через границы ячейки
308 286
 
309 287
 (defmethod border-flux ((grid grid2d-cv) (indx+ cell-index) (tau number) (border (eql :left)))
310 288
   (let* ((indx- (index (1- (x-index indx+))
311 289
 		       (y-index indx+)))
312
-	 (u+ (left-shot grid indx+))
313
-	 (u- (right-shot grid indx-))
290
+	 (u+ (primitive->conservative (left-shot (cell grid indx+))))
291
+	 (u- (primitive->conservative (right-shot (cell grid indx-))))
314 292
 	 (F- (flux u+))
315 293
 	 (F+ (flux u-))
316 294
 	 (a (pseudo-vel grid indx- 'right)))
... ...
@@ -323,100 +352,201 @@
323 323
 		     (y-index indx))))
324 324
     (border-flux grid indx tau :left)))
325 325
 
326
-;; солвер
327
-(defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number))
326
+(defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number) (method (eql :kurganov-tadmor-explicit)))
328 327
   (labels
329
-      ((predictor (grid tau)
328
+      ((predictor/corrector (grid tau)
330 329
 	 (let ((newgrid (clone grid)))
331 330
 	   (loop for j to (cells-max-y grid) do
332 331
 		(loop for i to (cells-max-x grid)
333 332
 		   for dF+ = (border-flux grid (index i j) tau :right)
334 333
 		   for dF- = (border-flux grid (index i j) tau :left) then dF+
335 334
 		   do
336
-		     (let* ((dF (cvar- dF+ dF-))
337
-			    (u (cell grid (index i j)))
338
-			    (hx (x-size u))
339
-			    (s (source u)))
340
-		       (labels (;; функция, для которой мы будем искать нули методом Ньютона
341
-				(eq-system (x1 x2 x3 x4)
342
-				  (let* ((x (make-instance 'conservative-variables
343
-							   :var-vector (vector x1 x2 x3 x4)))
344
-					 (result (cvar- (cvar+ (cvar/ (cvar- x u)
345
-								      tau)
346
-							       (cvar/ dF hx))
347
-							(cvar* 0.5 (cvar+ s (source x))))))
348
-				    (let-conservative result
349
-				      (values u1 u2 u3 u4))))
350
-				;; частные производеные этой функции (якобиан)
351
-				(eq-system-df (x1 x2 x3 x4)
352
-				  (let* ((D (+ (/ x2 x1)
353
-					       (/ x4 x3)))
354
-					 (F (* 0.5 *rl* D (signum D)))
355
-					 (F1 (/ (* F x2) (* x1 x1)))
356
-					 (F2 (/ F x1))
357
-					 (F3 (/ (* F x4) (* x3 x3)))
358
-					 (F4 (/ F x3))
359
-					 (TT (/ 1 tau)))
360
-				      (values TT 0.0 0.0 0.0 ;;
361
-					      F1 (- TT F2) (- F3) F4 ;;
362
-					      0.0 0.0 TT 0.0 ;;
363
-					      (- F1) F2 F3 (- TT F4))))
364
-				;; собирательная функция (необходима для GSLL)
365
-				(eq-system-fdf (x1 x2 x3 x4)
366
-				  (multiple-value-bind (v1 v2 v3 v4)
367
-				      (eq-system x1 x2 x3 x4)
368
-				    (multiple-value-bind (j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)
369
-					(eq-system-df x1 x2 x3 x4)
370
-				      (values v1 v2 v3 v4 j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)))))
371
-			 (let-conservative u
372
-			   (let ((max-iter *multi-root-iter*)
373
-				 (solver (make-multi-dimensional-root-solver-fdf
374
-					  +newton-mfdfsolver+
375
-					  (list #'eq-system #'eq-system-df #'eq-system-fdf)
376
-					  (grid:make-foreign-array
377
-					   'double-float :initial-contents (list u1 u2 u3 u4)))))
378
-			     (loop for iter from 0
379
-				while (if (< iter max-iter)
380
-					  (or (zerop iter)
381
-					      (not (multiroot-test-residual solver *multi-root-est-error*))))
382
-				do (iterate solver)
383
-				finally
384
-				  ;;(format t "solution: ~a~%" (solution solver)
385
-				  (setf (cell newgrid (index i j))
386
-					      (make-instance 'conservative-variables
387
-							     :var-vector (vector-double-float->vector
388
-									  (solution solver)))))))))))
389
-	   newgrid))
390
-       (corrector (grid predgrid tau)
391
-	 (let ((midgrid (clone grid))
392
-	       (newgrid (clone grid)))
393
-	   (loop for j to (cells-max-y grid) do
394
-		(loop for i to (cells-max-x grid) do
395
-		     (setf (cell midgrid (index i j))
396
-			   (cvar* 0.5 (cvar+ (cell grid (index i j))
397
-					     (cell predgrid (index i j)))))))
398
-	   (loop for j to (cells-max-y grid) do
399
-		(loop for i to (cells-max-x grid)
400
-		   for index = (index i j)
401
-		   for dF+ = (border-flux midgrid index tau :right)
402
-		   for dF- = (border-flux midgrid index tau :left) then dF+
403
-		   do
404
-		     (let* ((u (cell grid index))
405
-			    (S (cvar* 0.5 (cvar+ (source u)
406
-						 (source (cell predgrid index)))))
335
+		     (let* ((u (cell grid (index i j)))
336
+			    (S (source u))
407 337
 			    (hx (x-size u))
408 338
 			    (dF (cvar- dF+ dF-)))
409
-		       (setf (cell newgrid index)
339
+		       (setf (cell newgrid (index i j))
410 340
 			     (cvar- (cvar+ u (cvar* tau S))
411 341
 				    (cvar/ (cvar* tau dF)
412
-					   hx))))))
342
+					   hx)
343
+				    (cvar* tau (As*dU/dx grid (index i j))))))))
413 344
 	   newgrid)))
414
-    (let ((predgrid (predictor grid tau)))
415
-      (setf newgrid predgrid)
416
-      ;;(corrector grid predgrid tau)
417
-
418
-      (format t "PREDICTOR~%~%")
419
-      (print-grid newgrid :print-out-border t)
420
-      (newline)(newline)
421
-      )))
345
+    ;; считаем градиенты исходной сетки
346
+    (update-gradients grid)
347
+
348
+    (when (>= *debug-level* 10)
349
+      (format t "5: update-gradients~%")
350
+      (print-grid grid :print-out-border t)
351
+      (format t "----~%~%"))
352
+
353
+    ;; считаем приблизительные значения в (n+1)-й момент времени
354
+    (let ((predgrid (predictor/corrector grid tau)))
355
+
356
+      (when (>= *debug-level* 10)
357
+	(format t "6: predgrid-before~%")
358
+	(print-grid predgrid :print-out-border t)
359
+	(format t "----~%~%"))
360
+
361
+      ;; в качестве сетки, полученной на шаге-предикторе, берём полусумму старой и новой сеток
362
+      (loop for i to (cells-max-x grid) do
363
+	   (loop for j to (cells-max-y grid) do
364
+		(let* ((index (index i j))
365
+		       (old-vars (primitive->conservative (cell grid index)))
366
+		       (new-vars (primitive->conservative (cell predgrid index))))
367
+		  (setf (cell predgrid index)
368
+			(cvar* 0.5d0 (primitive->conservative 
369
+				      (cvar+ old-vars new-vars)))))))
370
+
371
+      (when (>= *debug-level* 10)
372
+	(format t "7: predgrid-after~%")
373
+	(print-grid predgrid :print-out-border t)
374
+	(format t "----~%~%"))
375
+      
376
+      ;; а градиенты оставляем теми же, что были на предикторе
377
+      (copy-gradients grid predgrid)
378
+
379
+      (when (>= *debug-level* 10)
380
+	(format t "8: copy-gradients~%")
381
+	(print-grid predgrid :print-out-border t)
382
+	(format t "----~%~%"))
383
+      
384
+      ;; повторяем тот же шаг, но с сеткой-предиктором и старыми коэффиэнтами - это и есть корректор
385
+      (setf newgrid (predictor/corrector predgrid tau))
386
+
387
+      (when (>= *debug-level* 10)
388
+	(format t "9: newgrid~%")
389
+	(print-grid newgrid :print-out-border t)
390
+	(format t "----~%~%"))
391
+      )
392
+
393
+    ;;newgrid ; в будущем планирую возвращать сетку, а не писваивать значения. Может зря.
394
+    ))
395
+
396
+;; солвер для неявной схемы курганова-тадмора
397
+;(defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number) (method (eql :kurganov-tadmor)))
398
+;  (labels
399
+;      ((predictor (grid tau)
400
+;	 (let ((newgrid (clone grid)))
401
+;	   (loop for j to (cells-max-y grid) do
402
+;		(loop for i to (cells-max-x grid)
403
+;		   for dF+ = (border-flux grid (index i j) tau :right)
404
+;		   for dF- = (border-flux grid (index i j) tau :left) then dF+
405
+;		   do
406
+;		     (let* ((dF (cvar- dF+ dF-))
407
+;			    (u (cell grid (index i j)))
408
+;			    (hx (x-size u))
409
+;			    (s (source u)))
410
+;		       (labels (;; функция, для которой мы будем искать нули методом Ньютона
411
+;				(eq-system (x1 x2 x3 x4)
412
+;				  (let* ((x (make-instance 'conservative-variables
413
+;							   :var-vector (vector x1 x2 x3 x4)))
414
+;					 (result (cvar- (cvar+ (cvar/ (cvar- x u)
415
+;								      tau)
416
+;							       (cvar/ dF hx))
417
+;							(cvar* 0.5 (cvar+ s (source x))))))
418
+;				    (let-conservative result
419
+;				      (values u1 u2 u3 u4))))
420
+;				;; частные производеные этой функции (якобиан)
421
+;				(eq-system-df (x1 x2 x3 x4)
422
+;				  (let* ((D (+ (/ x2 x1)
423
+;					       (/ x4 x3)))
424
+;					 (F (* 0.5 *rl* D (signum D)))
425
+;					 (F1 (/ (* F x2) (* x1 x1)))
426
+;					 (F2 (/ F x1))
427
+;					 (F3 (/ (* F x4) (* x3 x3)))
428
+;					 (F4 (/ F x3))
429
+;					 (TT (/ 1 tau)))
430
+;				      (values TT 0.0 0.0 0.0 ;;
431
+;					      F1 (- TT F2) (- F3) F4 ;;
432
+;					      0.0 0.0 TT 0.0 ;;
433
+;					      (- F1) F2 F3 (- TT F4))))
434
+;				;; собирательная функция (необходима для GSLL)
435
+;				(eq-system-fdf (x1 x2 x3 x4)
436
+;				  (multiple-value-bind (v1 v2 v3 v4)
437
+;				      (eq-system x1 x2 x3 x4)
438
+;				    (multiple-value-bind (j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)
439
+;					(eq-system-df x1 x2 x3 x4)
440
+;				      (values v1 v2 v3 v4 j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)))))
441
+;			 
442
+;			 (let-conservative u
443
+;			   (let ((max-iter *multi-root-iter*)
444
+;				 (solver (make-multi-dimensional-root-solver-fdf
445
+;					  +newton-mfdfsolver+
446
+;					  (list #'eq-system #'eq-system-df #'eq-system-fdf)
447
+;					  (grid:make-foreign-array
448
+;					   'double-float :initial-contents (list u1 u2 u3 u4)))))
449
+;			     (loop for iter from 0
450
+;				while (if (< iter max-iter)
451
+;					  (or (zerop iter)
452
+;					      (not (multiroot-test-residual solver *multi-root-est-error*))))
453
+;				do
454
+;				  (format t "iter ~a: ~a~%" iter (vector-double-float->vector (solution solver)))
455
+;				  (iterate solver)
456
+;				finally
457
+;				  ;;(format t "solution: ~a~%" (solution solver)
458
+;				  (setf (cell newgrid (index i j))
459
+;					      (make-instance 'conservative-variables
460
+;							     :var-vector (vector-double-float->vector
461
+;									  (solution solver)))))))))))
462
+;	   newgrid))
463
+;
464
+;       (corrector (grid predgrid tau)
465
+;	 (let ((midgrid (clone grid))
466
+;	       (newgrid (clone grid)))
467
+;	   (loop for j to (cells-max-y grid) do
468
+;		(loop for i to (cells-max-x grid) do
469
+;		     (setf (cell midgrid (index i j))
470
+;			   (cvar* 0.5 (cvar+ (cell grid (index i j))
471
+;					     (cell predgrid (index i j)))))))
472
+;	   (loop for j to (cells-max-y grid) do
473
+;		(loop for i to (cells-max-x grid)
474
+;		   for index = (index i j)
475
+;		   for dF+ = (border-flux midgrid index tau :right)
476
+;		   for dF- = (border-flux midgrid index tau :left) then dF+
477
+;		   do
478
+;		     (let* ((u (cell grid index))
479
+;			    (S (cvar* 0.5 (cvar+ (source u)
480
+;						 (source (cell predgrid index)))))
481
+;			    (hx (x-size u))
482
+;			    (dF (cvar- dF+ dF-)))
483
+;		       (setf (cell newgrid index)
484
+;			     (cvar- (cvar+ u (cvar* tau S))
485
+;				    (cvar/ (cvar* tau dF)
486
+;					   hx))))))
487
+;	   newgrid)))
488
+;    (let ((predgrid (predictor grid tau)))
489
+;      (setf newgrid predgrid)
490
+;      ;;(corrector grid predgrid tau)
491
+;
492
+;      (format t "PREDICTOR~%~%")
493
+;      (print-grid newgrid :print-out-border t)
494
+;      (newline)(newline)
495
+;      )))
496
+       
497
+
498
+				      
499
+			      
500
+
501
+;(defmacro list-of-vars (prefix argn &rest body)
502
+;  (let* ((fmt (format nil "~a~~d" (string-upcase (any->string prefix))))
503
+;	 (symbols (mapcar #'intern (loop for i from 1 to (eval argn)
504
+;				      collect (format nil fmt i)))))
505
+;    symbols))
506
+;
507
+;(defun inc-vector-elem (vec i dx)
508
+;  (let ((vec (copy-seq vec)))
509
+;    (incf (aref vec i) dx)
510
+;    vec))
511
+;
512
+;(defun make-jacobian (func &rest argv)
513
+;  (let ((tt 1d-5)
514
+;	(argc (1- (length argv))))
515
+;    (format t "1~%")
516
+;    (multiple-value-bind ((list-of-vars "f" argc) (apply func argv))
517
+;	(loop for i to argc do
518
+;	     (append (loop for j to argc do
519
+;			  (collect (/ (- (apply func (inc-vector-elem argv j tt)) (apply func argv))
520
+;				      tt))))))))
521
+;	 
522
+       
422 523
 
... ...
@@ -110,7 +110,7 @@
110 110
       )))
111 111
 
112 112
 
113
-(defun test-kt-constant (&key (cells-nx 1) (cells-ny 1) (output-dir #P"~/prog/clfpv/test-KT-constant/"))
113
+(defun test-kt-constant (&key (cells-nx 2) (cells-ny 1) (cu 1) (output-dir #P"~/prog/clfpv/test-KT-constant/"))
114 114
   "Тест на постоянном поле"
115 115
   (ensure-directories-exist output-dir)
116 116
   (let ((*DEFAULT-PATHNAME-DEFAULTS* output-dir))
... ...
@@ -121,6 +121,12 @@
121 121
 			       :right-boundary-condition 'zeroe
122 122
 			       :bottom-boundary-condition 'zeroe
123 123
 			       :top-boundary-condition 'zeroe)))
124
+
125
+      (when (>= *debug-level* 10)
126
+	(format t "1: ~%")
127
+	(print-grid grid :print-out-border t)
128
+	(format t "----~%~%"))
129
+
124 130
       (init-grid grid
125 131
 		 :left-border-coord 0.0d0
126 132
 		 :right-border-coord 1.0d0
... ...
@@ -128,18 +134,53 @@
128 128
 		 :bottom-border-line (constant-function 0.0d0)
129 129
 		 :thickening-coef-x 1
130 130
 		 :thickening-coef-y 1)
131
+
132
+      (when (>= *debug-level* 10)
133
+	(format t "2: ~%")
134
+	(print-grid grid :print-out-border t)
135
+	(format t "----~%~%"))
136
+
131 137
       (set-constant-field-2d grid
132 138
 			     (make-instance 'conservative-variables
133 139
 					    :phi 0.999999d0
134 140
 					    :P 100000.0d0
135 141
 					    :ug 25.0d0
136 142
 					    :ul 25.0d0))
143
+
144
+      (when (>= *debug-level* 10)
145
+	(format t "3: ~%")
146
+	(print-grid grid :print-out-border t)
147
+	(format t "----~%~%"))
148
+
137 149
       (update-boundary-cells grid)
138
-      (let ((newgrid (march grid :finish-step 1)))
139
-	newgrid
140
-	))))
141 150
 
151
+      (when (>= *debug-level* 10)
152
+	(format t "4.1: ~%")
153
+	(print-grid grid :print-out-border t)
154
+	(format t "----~%~%"))
155
+      
156
+      ;;(update-gradients grid)
157
+
158
+      (when (>= *debug-level* 10)
159
+	(format t "4.2: ~%")
160
+	(print-grid grid :print-out-border t)
161
+	(format t "----~%~%"))
142 162
 
143 163
 
164
+      (when (>= *debug-level* 3)
165
+	(format t "start-grid: ~%")
166
+	(print-grid grid :print-out-border t)
167
+	(format t "----~%~%"))
168
+      
169
+      (let ((newgrid (march grid
170
+			    :finish-step 10
171
+			    :cu cu
172
+			    )))
144 173
 
145
-    
174
+	(when (>= *debug-level* 3)
175
+	  (format t "final-grid: ~%")
176
+	  (print-grid newgrid :print-out-border t)
177
+	  (format t "----~%~%"))
178
+	
179
+	newgrid
180
+	))))
... ...
@@ -161,3 +161,6 @@
161 161
 
162 162
 (defmethod minmod ((a conservative-variables) (b conservative-variables))
163 163
   (conservative-variables-operator (list a b) #'minmod))
164
+
165
+(defmethod minmod ((a primitive-variables) (b primitive-variables))
166
+  (primitive-variables-operator (list a b) #'minmod))
... ...
@@ -121,5 +121,10 @@
121 121
 (defun pvar-signum (pv)
122 122
   (primitive-variables-operator (list pv) #'signum))
123 123
 
124
+;;
124 125
 
125
-
126
+(defmethod make-dummy-variables ((type (eql 'primitive-variables)))
127
+  (make-instance 'primitive-variables
128
+		 :var-vector (coerce (loop repeat *number-of-primitive-variables*
129
+					collect 'undef)
130
+				     'vector)))