Browse code

Добавил дебаги. Заменил minmod-ы. Добился корректной работы test-kt-split-1

Dmitrii Kashin authored on 25/05/2015 17:37:09
Showing 8 changed files
... ...
@@ -26,7 +26,24 @@
26 26
 (defparameter *multi-root-est-error* 1.0d-8)
27 27
 
28 28
 (defparameter *constant-timestep* .5d-8)
29
+(defparameter *timestep-multiplier* 1.0d-4)
29 30
 
30 31
 (defparameter *source-timestep-multiplier* 0.2d0)
31 32
 
32
-(defparameter *debug-level* 5)
33
+(defparameter *minmod-coef* 1.5) ;; must be 1<=x<=2
34
+
35
+(defparameter *debug-level* 0)
36
+
37
+(defparameter *debug-shots* nil)
38
+
39
+(defparameter *debug-newton* nil)
40
+(defparameter *debug-grids-in-make-step* nil)
41
+(defparameter *debug-radius* nil)
42
+
43
+(defparameter *debug-messages* nil)
44
+(defparameter *debug-predictor/corrector* nil)
45
+(defparameter *debug-pseudo-vel* nil)
46
+
47
+
48
+(defparameter *enable-breaks* nil)
49
+
... ...
@@ -33,3 +33,19 @@
33 33
 		    (if (< i mid-x-index)
34 34
 			left-vars
35 35
 			right-vars))))))
36
+
37
+(defun perturbate-field-x (grid perturbate)
38
+  "Эта функция применяет к каждой ячейке сетки `grid' с установленным
39
+полем функцию `perturbate', которая вызывается с индексом по x текущей
40
+ячейки. Параметры этой ячейки в сетке домножаются на полученный
41
+коэффициент."
42
+  (let* ((cells-amount-x (cells-amount-x grid))
43
+	 (cells-amount-y (cells-amount-y grid))
44
+	 (cells-max-x (1- cells-amount-x))
45
+	 (cells-max-y (1- cells-amount-y)))
46
+    (loop for i from 0 to cells-max-x do
47
+	 (loop for j from 0 to cells-max-y do
48
+	      (let* ((pv (conservative->primitive (cell grid (index i j))))
49
+		     (pv (pvar* pv (funcall perturbate i))))
50
+		(setf (var-vector (cell grid (index i j)))
51
+		      (var-vector (primitive->conservative pv))))))))
... ...
@@ -78,6 +78,15 @@
78 78
 
79 79
 ;; вспомогательные функции
80 80
 
81
+(defun enbreak (text)
82
+  "Остановиться и запустить дебаггер, если выставлен в истину флаг
83
+`*enable-breaks*'"
84
+  (when *enable-breaks* (break text)))
85
+
86
+(defun dbgmsg (text)
87
+  (when *debug-messages*
88
+    (format t text)))
89
+
81 90
 (defun vector-double-float->vector (agrid)
82 91
   "В документации сказано (см. раздел 11 системы antik), что объекты
83 92
 типа `vector-double-float' можно преобразовать к типу `vector'
... ...
@@ -299,9 +308,5 @@
299 299
   (:documentation "Устанавливает в сетке `grid' заграничные ячейки в
300 300
   соответствии с её граничными условиями"))
301 301
 
302
-(defgeneric minmod (a b)
302
+(defgeneric minmod (type &rest params)
303 303
   (:documentation "Вычисляет минимум по модулю"))
304
-
305
-
306
-
307
-
... ...
@@ -13,7 +13,7 @@
13 13
   (unless (or finish-step finish-time)
14 14
     (error "По крайней мере одно условие остановки должно быть указано"))
15 15
   (when finish-step
16
-    (format nil "Остановка после ~d итераций.~%" finish-step))
16
+    (format t "Остановка после ~d итераций.~%" finish-step))
17 17
   (when finish-time
18 18
     (format t "Остановка после ~a секунд.~%" finish-time))
19 19
 
... ...
@@ -23,7 +23,7 @@
23 23
     ;; сначала экспортируем начальное поле
24 24
     (export-grid grid :print-out-border print-out-border)
25 25
     ;; если finish-step установлен в 0, то больше ничего не делаем
26
-    (when (<= finish-step 0)
26
+    (when (and finish-step (<= finish-step 0))
27 27
       (setf stop-flag t))
28 28
     (loop until stop-flag do
29 29
 	 (update-gradients grid)
... ...
@@ -57,7 +57,8 @@
57 57
 		      (zerop (mod (1+ iterations) output-every-step)))
58 58
 	     (setf print-flag t))
59 59
 	   ;; проверяем, является ли следующий шаг последним
60
-	   (when (>= (1+ iterations) finish-step)
60
+	   (when (and finish-step
61
+		      (>= (1+ iterations) finish-step))
61 62
 	     (setf print-flag t)
62 63
 	     (setf stop-flag t))
63 64
 	   
... ...
@@ -23,3 +23,9 @@
23 23
 (defun export-grid&nodes (grid)
24 24
   (export-grid grid :print-out-border t)
25 25
   (export-nodes grid))
26
+
27
+(defun gridbg (text grid)
28
+  (when *debug-grids-in-make-step*
29
+    (format t "~a: ~%" text)
30
+    (print-grid grid :print-out-border nil)
31
+    (format t "----~%~%")))
... ...
@@ -78,17 +78,31 @@
78 78
 		    (values (frolov-equation x)
79 79
 			    (frolov-equation-derivative x)))
80 80
 		  (find-root (&optional initial (method +newton-fdfsolver+))
81
+		    (when *debug-newton*
82
+		      (format t "~%cell:~%")
83
+		      (print-content t cell)
84
+		      (newline))
81 85
 		    (let* ((max-iter *newton-max-iter*)
82 86
 			   (solver (make-one-dimensional-root-solver-fdf
83 87
 				    method #'frolov-equation #'frolov-equation-derivative #'frolov-equation-fdf initial)))
84 88
 		      (loop for iter from 0
85 89
 			 for oldroot = initial then root
86 90
 			 for root = (progn (iterate solver) (solution solver))
91
+			 for olderror = 0 then error
92
+			 for error = (- root oldroot)
87 93
 			 while (if (< iter max-iter)
88 94
 				   (not (root-test-delta root oldroot *newton-est-error* 0.0d0))
89 95
 				   (error "Метод Ньютона превысил допустимое число итераций."))
96
+			 do
97
+			   (when *debug-newton*
98
+			     (when (zerop iter)
99
+			       (format t "NEWTON: initial-root ~d~%" oldroot))
100
+			     (format t "NEWTON: iter ~d, root ~d, error ~d~%"
101
+				     (1+ iter) root error))
102
+			   (when (minusp (* olderror error))
103
+			     (error "Нарушена монотонность метода Ньютона. Возможно, получилось отрицательное давление?"))
90 104
 			 finally (return root)))))
91
-	   (let* ((suggestion (* 5 cg))
105
+	   (let* ((suggestion (+ (* 5 cg) 20000)) ;; 20000 - взято от балды
92 106
 		  (left-init (- (min ul (- ug cg) (+ ug cg))
93 107
 				suggestion))
94 108
 		  (right-init (+ (max ul (- ug cg) (+ ug cg))
... ...
@@ -163,9 +177,13 @@
163 163
 ;;(eigen-values-by-cell *eigen-cell* 'S) => -3e-7
164 164
 ;;(eigen-values-by-cell *eigen-cell* 'Afull)
165 165
 
166
-;; спектральный радиус
167 166
 (defun radius (cell type)
168
-  (apply #'max (mapcar #'abs (eigen-values-by-cell cell type))))
167
+  "Спектральный радиус в методе Курганова-Тадмора"
168
+  (let ((radius (apply #'max (mapcar #'abs (eigen-values-by-cell cell type)))))
169
+    (when *debug-radius*
170
+      (format t "RADIUS: ~a~%" radius))
171
+    radius))
172
+      
169 173
 
170 174
 ;; потоковый член в ячейке.
171 175
 (defmethod flux ((vars conservative-variables))
... ...
@@ -244,12 +262,19 @@
244 244
 		  (xm (x-coord mid-cell))
245 245
 		  (xl (unless outborder-left (x-coord left-cell)))
246 246
 		  (xr (unless outborder-right (x-coord right-cell))))
247
-	     (minmod (pvar/ (pvar- (conservative->primitive right-cell)
247
+	     (minmod 'primitive-variables
248
+		     (pvar* *minmod-coef*
249
+			    (pvar/ (pvar- (conservative->primitive right-cell)
250
+					  (conservative->primitive mid-cell))
251
+				   (- xr xm)))
252
+		     (pvar* *minmod-coef*
253
+			    (pvar/ (pvar- (conservative->primitive mid-cell)
254
+					  (conservative->primitive left-cell))
255
+				   (- xm xl)))
256
+		     (pvar/ (pvar- (conservative->primitive right-cell)
248 257
 				   (conservative->primitive mid-cell))
249
-			    (- xr xm))
250
-		     (pvar/ (pvar- (conservative->primitive mid-cell)
251
-				   (conservative->primitive left-cell))
252
-			    (- xm xl))))))))
258
+			    (* 2 (- xr xl)))
259
+		     ))))))
253 260
 
254 261
 (defmethod update-gradients ((grid grid2d-cv))
255 262
   (loop for i from (min-x-index/border grid) to (max-x-index/border grid) do
... ...
@@ -300,11 +325,20 @@
300 300
   (let ((pvars (conservative->primitive cell))
301 301
 	(gx (gradients cell))
302 302
 	(hx (x-size cell)))
303
-    (primitive->conservative
304
-     (case direction
305
-       ((left) (pvar- pvars (pvar* 0.5 hx gx)))
306
-       ((right) (pvar+ pvars (pvar* 0.5 hx gx)))
307
-       (otherwise (error 'shot "Неизвестное направление для стрельбы"))))))
303
+    (let ((shot (primitive->conservative
304
+		 (case direction
305
+		   ((left) (pvar- pvars (pvar* 0.5 hx gx)))
306
+		   ((right) (pvar+ pvars (pvar* 0.5 hx gx)))
307
+		   (otherwise (error 'shot "Неизвестное направление для стрельбы"))))))
308
+      (when *debug-shots*
309
+	(format t "~%SHOT: ~a~%" direction)
310
+	(format t "SHOT: from-cell:~%")
311
+	(print-content t cell)
312
+	(format t "~%SHOT: result:~%")
313
+	(print-content t shot)
314
+	(newline)
315
+	)
316
+      shot)))
308 317
 
309 318
 (defun left-shot (cell)
310 319
   (shot cell 'left))
... ...
@@ -323,15 +357,20 @@
323 323
 	 (mid-cell (cell grid indx))
324 324
 	 (left-cell (unless outborder-left (cell grid l-indx)))
325 325
 	 (right-cell (unless outborder-right (cell grid r-indx))))
326
-    (case direction
327
-      ((right) (max (radius (primitive->conservative (right-shot mid-cell)) 'A)
328
-		    (radius (primitive->conservative (left-shot right-cell)) 'A)
329
-		    (radius mid-cell 'Afull)
330
-		    (radius right-cell 'Afull)))
331
-      ((left) (max (radius (primitive->conservative (right-shot left-cell)) 'A)
332
-		   (radius (primitive->conservative (left-shot mid-cell)) 'A)
333
-		   (radius left-cell 'Afull)
334
-		   (radius right-cell 'Afull))))))
326
+    (when *debug-newton*
327
+      (format t "~%NEWTON: Считаем псевдоскорость в ячейке (~d ~d)~%" i j))
328
+    (let ((vel (case direction
329
+		 ((right) (max (radius (primitive->conservative (right-shot mid-cell)) 'A)
330
+			       (radius (primitive->conservative (left-shot right-cell)) 'A)
331
+			       (radius mid-cell 'Afull)
332
+			       (radius right-cell 'Afull)
333
+			       ))
334
+		 ((left) (max (radius (primitive->conservative (right-shot left-cell)) 'A)
335
+			      (radius (primitive->conservative (left-shot mid-cell)) 'A)
336
+			      (radius left-cell 'Afull)
337
+			      (radius right-cell 'Afull))))))
338
+      (when *debug-pseudo-vel* (format t "~%PSEUDO-VEL: ~d~%" vel))
339
+      vel)))
335 340
 
336 341
 ;; потоки через границы ячейки
337 342
 
... ...
@@ -344,8 +383,8 @@
344 344
 	 (F+ (flux u-))
345 345
 	 (a (pseudo-vel grid indx- 'right)))
346 346
     (cvar* 0.5
347
-	   (cvar- (cvar+ F+ F-))
348
-	   (cvar* a (cvar- u+ u-)))))
347
+	   (cvar- (cvar+ F+ F-)
348
+		  (cvar* a (cvar- u+ u-))))))
349 349
 
350 350
 (defmethod border-flux ((grid grid2d-cv) (indx cell-index) (tau number) (border (eql :right)))
351 351
   (let ((indx (index (1+ (x-index indx))
... ...
@@ -354,39 +393,48 @@
354 354
 
355 355
 (defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number) (method (eql :kurganov-tadmor-explicit)))
356 356
   (labels
357
-      ((predictor/corrector (grid tau)
358
-	 (let ((newgrid (clone grid)))
359
-	   (loop for j to (cells-max-y grid) do
360
-		(loop for i to (cells-max-x grid)
361
-		   for dF+ = (border-flux grid (index i j) tau :right)
362
-		   for dF- = (border-flux grid (index i j) tau :left) then dF+
363
-		   do
364
-		     (let* ((u (cell grid (index i j)))
365
-			    (S (source u))
366
-			    (hx (x-size u))
367
-			    (dF (cvar- dF+ dF-)))
368
-		       (setf (cell newgrid (index i j))
369
-			     (cvar- (cvar+ u (cvar* tau S))
370
-				    (cvar/ (cvar* tau dF)
371
-					   hx)
372
-				    (cvar* tau (As*dU/dx grid (index i j))))))))
373
-	   newgrid)))
357
+      ((predictor/corrector (grid newgrid tau)
358
+	 (loop for j to (cells-max-y grid) do
359
+	      (loop for i to (cells-max-x grid)
360
+		 for dF- = (border-flux grid (index i j) tau :left) then dF+
361
+		 for dF+ = (border-flux grid (index i j) tau :right)
362
+		 do
363
+		   (let* ((u (cell grid (index i j)))
364
+			  (S (source u))
365
+			  (hx (x-size u))
366
+			  (dF (cvar- dF+ dF-))
367
+			  (P*Dphi (As*dU/dx grid (index i j)))
368
+			  (nu (cvar- (cvar+ u (cvar* tau S))
369
+				     (cvar/ (cvar* tau dF)
370
+					    hx)
371
+				     (cvar* tau P*Dphi))))
372
+		     (when *debug-predictor/corrector*
373
+		       (format t "~%~@{~{PRED-CORR: ~a = ~d~%~}~}"
374
+			       `((i j) (,i ,j))
375
+			       `(oldcell-cv ,(print-conservative nil u :with-headers nil))
376
+			       `(oldcell-pv ,(print-content nil u :with-headers nil))
377
+			       `(tau ,tau)
378
+			       `(dF+ ,(print-conservative nil dF+ :with-headers nil))
379
+			       `(dF- ,(print-conservative nil dF- :with-headers nil))
380
+			       `(dF ,(print-conservative nil dF :with-headers nil))
381
+			       `(S ,(print-conservative nil S :with-headers nil))
382
+			       `(P*Dphi ,(print-conservative nil P*Dphi :with-headers nil))
383
+			       `(newcell-cv ,(print-conservative nil nu :with-headers nil))
384
+			       `(newcell-pv ,(print-content nil nu :with-headers nil)))
385
+		       (enbreak "Шаг предиктора/корректора"))
386
+		     (setf (cell newgrid (index i j)) nu)
387
+		     ))
388
+	      (gridbg "X: predictor/corrector~%" newgrid))))
374 389
     ;; считаем градиенты исходной сетки
375 390
     (update-gradients grid)
376
-
377
-    (when (>= *debug-level* 10)
378
-      (format t "5: update-gradients~%")
379
-      (print-grid grid :print-out-border t)
380
-      (format t "----~%~%"))
381
-
391
+    (gridbg "5: update-gradients~%" grid)
382 392
     ;; считаем приблизительные значения в (n+1)-й момент времени
383
-    (let ((predgrid (predictor/corrector grid tau)))
384
-
385
-      (when (>= *debug-level* 10)
386
-	(format t "6: predgrid-before~%")
387
-	(print-grid predgrid :print-out-border t)
388
-	(format t "----~%~%"))
389
-
393
+    (dbgmsg "------------------------------------------------------------~%")
394
+    (dbgmsg "Делаем предиктор~%")
395
+    (when (or *debug-shots* *debug-predictor/corrector*) (print-grid grid :print-out-border nil))
396
+    (let ((predgrid (clone grid)))
397
+      (predictor/corrector grid predgrid tau)
398
+      (gridbg "6: predgrid-before~%" predgrid)
390 399
       ;; в качестве сетки, полученной на шаге-предикторе, берём полусумму старой и новой сеток
391 400
       (loop for i to (cells-max-x grid) do
392 401
 	   (loop for j to (cells-max-y grid) do
... ...
@@ -396,29 +444,20 @@
396 396
 		  (setf (cell predgrid index)
397 397
 			(cvar* 0.5d0 (primitive->conservative 
398 398
 				      (cvar+ old-vars new-vars)))))))
399
-
400
-      (when (>= *debug-level* 10)
401
-	(format t "7: predgrid-after~%")
402
-	(print-grid predgrid :print-out-border t)
403
-	(format t "----~%~%"))
404
-      
399
+      (gridbg "7: predgrid-after~%" predgrid)
405 400
       ;; а градиенты оставляем теми же, что были на предикторе
406 401
       (copy-gradients grid predgrid)
407
-
408
-      (when (>= *debug-level* 10)
409
-	(format t "8: copy-gradients~%")
410
-	(print-grid predgrid :print-out-border t)
411
-	(format t "----~%~%"))
412
-      
402
+      (gridbg "8: copy-gradients~%" predgrid)
403
+      ;; надо удовлетворять граничным условиям
404
+      (update-boundary-cells predgrid)
405
+      (gridbg "8.1: update-boundary-cells~%" predgrid)
413 406
       ;; повторяем тот же шаг, но с сеткой-предиктором и старыми коэффиэнтами - это и есть корректор
414
-      (setf newgrid (predictor/corrector predgrid tau))
415
-
416
-      (when (>= *debug-level* 10)
417
-	(format t "9: newgrid~%")
418
-	(print-grid newgrid :print-out-border t)
419
-	(format t "----~%~%"))
407
+      (dbgmsg "------------------------------------------------------------~%")
408
+      (dbgmsg "Делаем корректор~%")
409
+      (when (or *debug-shots* *debug-predictor/corrector*) (print-grid predgrid :print-out-border nil))
410
+      (predictor/corrector predgrid newgrid tau)
411
+      (gridbg "9: newgrid~%" newgrid)
420 412
       )
421
-
422 413
     ;;newgrid ; в будущем планирую возвращать сетку, а не писваивать значения. Может зря.
423 414
     ))
424 415
 
... ...
@@ -1,12 +1,6 @@
1 1
 (in-package :clf)
2 2
 
3
-(defun kt-dbg (level text grid)
4
-  (when (>= *debug-level* level)
5
-    (format t "~a: ~%" text)
6
-    (print-grid grid :print-out-border nil)
7
-    (format t "----~%~%")))
8
-
9
-(defun test-kt-constant (&key (cells-nx 10) (cells-ny 1) (cu 1) (output-dir #P"~/prog/clfpv/test-KT-constant/"))
3
+(defun test-kt-constant (&key (cells-nx 4) (cells-ny 1) (cu 1) (output-dir #P"~/prog/clfpv/test-kt-constant/"))
10 4
   "Тест на постоянном поле. Он нужен в основном для того, чтобы
11 5
 проверить консервативность схемы, а также для дебага программы."
12 6
   (ensure-directories-exist output-dir)
... ...
@@ -18,7 +12,7 @@
18 18
 			       :right-boundary-condition 'zeroe
19 19
 			       :bottom-boundary-condition 'zeroe
20 20
 			       :top-boundary-condition 'zeroe)))
21
-      (kt-dbg 10 "1" grid)
21
+      (gridbg "1" grid)
22 22
       (init-grid grid
23 23
 		 :left-border-coord 0.0d0
24 24
 		 :right-border-coord 1.0d0
... ...
@@ -26,28 +20,30 @@
26 26
 		 :bottom-border-line (constant-function 0.0d0)
27 27
 		 :thickening-coef-x 1
28 28
 		 :thickening-coef-y 1)
29
-      (kt-dbg 10 "2" grid)
29
+      (gridbg "2" grid)
30 30
       (set-constant-field-2d grid
31 31
 			     (make-instance 'conservative-variables
32 32
 					    :phi 0.999999d0
33 33
 					    :P 100000.0d0
34 34
 					    :ug 25.0d0
35 35
 					    :ul 25.0d0))
36
-      (kt-dbg 10 "3" grid)
36
+      (gridbg "3" grid)
37 37
       (update-boundary-cells grid)
38
-      (kt-dbg 10 "4" grid)
39
-      (kt-dbg 3 "start-grid" grid)
38
+      (gridbg "4" grid)
39
+      (gridbg "start-grid" grid)
40 40
       (let ((newgrid (march grid
41 41
 			    :finish-step 10
42
+			    :output-every-step 1
42 43
 			    :cu cu
43 44
 			    )))
44
-	(kt-dbg 3 "final-grid" newgrid)
45
+	(gridbg "final-grid" newgrid)
45 46
 	newgrid
46 47
 	))))
47 48
 
48
-(defun test-kt-split-1 (&key (cells-nx 100) (cells-ny 1) (cu 1) (output-dir #P"~/prog/clfpv/test-kt-split-1/"))
49
-  "Тест на скачке. В данном тесте значение phi близко к 1, а это
50
-значит, что решение должно быть почти таким же, как в газе."
49
+(defun test-kt-split-1 (&key (cells-nx 100) (cells-ny 1) (cu 0.001) (output-dir #P"~/prog/clfpv/test-kt-split-1/"))
50
+  "Тест на стоячем скачке уплотнения. В данном тесте значение phi
51
+близко к 1, а это значит, что решение должно быть почти таким же, как
52
+в газе."
51 53
   (ensure-directories-exist output-dir)
52 54
   (let ((*DEFAULT-PATHNAME-DEFAULTS* output-dir))
53 55
     (let ((grid (make-instance 'grid2d-cv
... ...
@@ -57,7 +53,7 @@
57 57
 			       :right-boundary-condition 'zeroe
58 58
 			       :bottom-boundary-condition 'zeroe
59 59
 			       :top-boundary-condition 'zeroe)))
60
-      (kt-dbg 10 "1 (instance)" grid)
60
+      (gridbg "1 (instance)" grid)
61 61
       (init-grid grid
62 62
 		 :left-border-coord 0.0d0
63 63
 		 :right-border-coord 1.0d0
... ...
@@ -65,7 +61,7 @@
65 65
 		 :bottom-border-line (constant-function 0.0d0)
66 66
 		 :thickening-coef-x 1
67 67
 		 :thickening-coef-y 1)
68
-      (kt-dbg 10 "2 (init-grid)" grid)
68
+      (gridbg "2 (init-grid)" grid)
69 69
       (let* ((left-vars (make-instance 'primitive-variables
70 70
 				       :phi 0.999999
71 71
 				       :rg 1.0
... ...
@@ -84,16 +80,198 @@
84 84
 	(set-splitted-field-2d-x grid
85 85
 				 :left-vars left-vars
86 86
 				 :right-vars right-vars))
87
-      (kt-dbg 10 "3 (set-field)" grid)
87
+      (gridbg "3 (set-field)" grid)
88
+      (update-boundary-cells grid)
89
+      (gridbg "4 (boundaries)" grid)
90
+      (gridbg "start-grid" grid)
91
+      (let ((newgrid (march grid
92
+			    :finish-time 0.05
93
+			    :output-every-time 0.0001
94
+			    :cu cu
95
+			    )))
96
+	(gridbg "final-grid" newgrid)
97
+	newgrid
98
+	)
99
+      )))
100
+
101
+(defun test-kt-split-1v2 (&key (cells-nx 100) (cells-ny 1) (cu 0.001) (output-dir #P"~/prog/clfpv/test-kt-split-1v2/"))
102
+  "Тест на стоячем скачке уплотнения. В данном тесте значение phi
103
+близко к 1, а это значит, что решение должно быть почти таким же, как
104
+в газе. Этот тест отличается от предыдущего тем, что во всей сетке
105
+скорости увеличены, чтобы скачок двигался. Таким образом можно
106
+наблюдать размазывание."
107
+  (ensure-directories-exist output-dir)
108
+  (let ((*DEFAULT-PATHNAME-DEFAULTS* output-dir))
109
+    (let ((grid (make-instance 'grid2d-cv
110
+			       :cells-amount-x cells-nx
111
+			       :cells-amount-y cells-ny
112
+			       :left-boundary-condition 'zeroe
113
+			       :right-boundary-condition 'zeroe
114
+			       :bottom-boundary-condition 'zeroe
115
+			       :top-boundary-condition 'zeroe)))
116
+      (gridbg "1 (instance)" grid)
117
+      (init-grid grid
118
+		 :left-border-coord 0.0d0
119
+		 :right-border-coord 1.0d0
120
+		 :top-border-line (constant-function 1.0d0)
121
+		 :bottom-border-line (constant-function 0.0d0)
122
+		 :thickening-coef-x 1
123
+		 :thickening-coef-y 1)
124
+      (gridbg "2 (init-grid)" grid)
125
+      (let* ((left-vars (make-instance 'primitive-variables
126
+				       :phi 0.999999
127
+				       :rg 1.0
128
+				       :ug 400.0
129
+				       :ul 400.0))
130
+	     (right-vars (make-instance 'primitive-variables
131
+					:phi (phi left-vars)
132
+					:ug (/ (* *sos* *sos*)
133
+					       (* *kappa* (ug left-vars)))
134
+					:ul (/ (* *sos* *sos*)
135
+					       (* *kappa* (ul left-vars)))
136
+					:rg (/ (* (rg left-vars) (ug left-vars))
137
+					       ;; на ug справа
138
+					       (/ (* *sos* *sos*)
139
+						  (* *kappa* (ug left-vars)))))))
140
+	(set-splitted-field-2d-x grid
141
+				 :left-vars left-vars
142
+				 :right-vars right-vars))
143
+      (loop for i from (min-x-index grid) to (max-x-index grid) do
144
+	   (loop for j from (min-y-index grid) to (min-y-index grid) do
145
+		(let ((pv (conservative->primitive (cell grid (index i j)))))
146
+		  (incf (ug pv) 100)
147
+		  (incf (ul pv) 100)
148
+		  (setf (cell grid (index i j)) (primitive->conservative pv)))))
149
+      (format t "~%INCREMENT DEBUG:~%")
150
+      (print-content t (cell grid (index 0 0)))
151
+      
152
+      (gridbg "3 (set-field)" grid)
153
+      (update-boundary-cells grid)
154
+      (gridbg "4 (boundaries)" grid)
155
+      (gridbg "start-grid" grid)
156
+      (let ((newgrid (march grid
157
+			    :finish-step 1
158
+			    :output-every-step 1
159
+			    :cu cu
160
+			    )))
161
+	(gridbg "final-grid" newgrid)
162
+	newgrid
163
+	)
164
+      )))
165
+
166
+(defun test-kt-split-2 (&key (cells-nx 100) (cells-ny 1) (cu 1) (output-dir #P"~/prog/clfpv/test-kt-split-2/"))
167
+  "Тест на стоячем скачке уплотнения при phi = 0.9. Этот тест нужен,
168
+чтобы проверить, как ведёт себя метод при расчёте не чисто газового
169
+течения."
170
+  (ensure-directories-exist output-dir)
171
+  (let ((*DEFAULT-PATHNAME-DEFAULTS* output-dir))
172
+    (let ((grid (make-instance 'grid2d-cv
173
+			       :cells-amount-x cells-nx
174
+			       :cells-amount-y cells-ny
175
+			       :left-boundary-condition 'zeroe
176
+			       :right-boundary-condition 'zeroe
177
+			       :bottom-boundary-condition 'zeroe
178
+			       :top-boundary-condition 'zeroe)))
179
+      (gridbg "1 (instance)" grid)
180
+      (init-grid grid
181
+		 :left-border-coord 0.0d0
182
+		 :right-border-coord 1.0d0
183
+		 :top-border-line (constant-function 1.0d0)
184
+		 :bottom-border-line (constant-function 0.0d0)
185
+		 :thickening-coef-x 1
186
+		 :thickening-coef-y 1)
187
+      (gridbg "2 (init-grid)" grid)
188
+
189
+      (let* ((left (make-instance 'primitive-variables
190
+				  :phi 0.9
191
+				  :rg 1.0
192
+				  :ug 1000.0
193
+				  :ul 1000.0))
194
+	     (ug-right (ug left))
195
+	     (ul-right (/ (+ (* (rg left) *Ro* *T*)
196
+			     (* (- 1 (phi left)) *rl* (ul left) (ul left)))
197
+			  (* *rl* (ul left))))
198
+	     (phi-right (- 1
199
+			   (/ (* (- 1 (phi left)) (ul left))
200
+			      ul-right)))
201
+	     (rg-right (/ (* (rg left) (phi left))
202
+			  phi-right))
203
+	     (right
204
+	      (make-instance 'primitive-variables
205
+			     :ug ug-right
206
+			     :ul ul-right
207
+			     :phi phi-right
208
+			     :rg rg-right)))
209
+	(set-splitted-field-2d-x grid
210
+				 :left-vars left
211
+				 :right-vars right))
212
+      
213
+      (gridbg "3 (set-field)" grid)
214
+      (update-boundary-cells grid)
215
+      (gridbg "4 (boundaries)" grid)
216
+      (gridbg "start-grid" grid)
217
+      (let ((newgrid (march grid
218
+			    :finish-step 1000
219
+			    :output-every-step 50
220
+			    :cu cu
221
+			    )))
222
+	(gridbg "final-grid" newgrid)
223
+	newgrid
224
+	)
225
+      )))
226
+
227
+
228
+
229
+(defun test-kt-frolov-b (&key (cells-nx 100) (cells-ny 1) (cu 1) (output-dir #P"~/prog/clfpv/test-kt-frolov-b/"))
230
+  "Тест из статьи Фролова с источниковыми членами. В статье обозначен как Б."
231
+  (ensure-directories-exist output-dir)
232
+  (let ((*DEFAULT-PATHNAME-DEFAULTS* output-dir))
233
+    (let ((grid (make-instance 'grid2d-cv
234
+			       :cells-amount-x cells-nx
235
+			       :cells-amount-y cells-ny
236
+			       :left-boundary-condition 'zeroe
237
+			       :right-boundary-condition 'zeroe
238
+			       :bottom-boundary-condition 'zeroe
239
+			       :top-boundary-condition 'zeroe)))
240
+      (gridbg "1 (instance)" grid)
241
+      (init-grid grid
242
+		 :left-border-coord 0.0d0
243
+		 :right-border-coord 1.0d0
244
+		 :top-border-line (constant-function 1.0d0)
245
+		 :bottom-border-line (constant-function 0.0d0)
246
+		 :thickening-coef-x 1
247
+		 :thickening-coef-y 1)
248
+      (gridbg "2 (init-grid)" grid)
249
+
250
+      (let ((Pin (* 20 100000.0d0))
251
+	    (Pout (* 10 100000.0d0))
252
+	    (phi 0.01d0)
253
+	    (ug 25.0d0)
254
+	    (ul 30.0d0))
255
+	(let ((left (make-instance 'primitive-variables
256
+				   :phi phi
257
+				   :P Pin
258
+				   :ug ug
259
+				   :ul ul))
260
+	      (right (make-instance 'primitive-variables
261
+				    :phi phi
262
+				    :P Pout
263
+				    :ug ug
264
+				    :ul ul)))
265
+	  (set-monotonous-field-2d-x grid
266
+				     :left-vars left
267
+				     :right-vars right)))
268
+      
269
+      (gridbg "3 (set-field)" grid)
88 270
       (update-boundary-cells grid)
89
-      (kt-dbg 10 "4 (boundaries)" grid)
90
-      (kt-dbg 3 "start-grid" grid)
271
+      (gridbg "4 (boundaries)" grid)
272
+      (gridbg "start-grid" grid)
91 273
       (let ((newgrid (march grid
92 274
 			    :finish-step 1000
93 275
 			    :output-every-step 50
94 276
 			    :cu cu
95 277
 			    )))
96
-	(kt-dbg 3 "final-grid" newgrid)
278
+	(gridbg "final-grid" newgrid)
97 279
 	newgrid
98 280
 	)
99 281
       )))
... ...
@@ -30,11 +30,10 @@
30 30
 
31 31
 ;; дополнительная возможность вывода в чистом виде, как они есть (для дебага)
32 32
 
33
-(defun print-conservative (stream cv)
34
-  (format stream "u1 u2 u3 u4")
35
-  (newline stream)
36
-  (apply #'format stream "~d ~d ~d ~d" (vector->list (var-vector cv)))
37
-  (newline stream))
33
+(defun print-conservative (stream cv &key (with-headers t))
34
+  (when with-headers
35
+    (format stream "u1 u2 u3 u4~%"))
36
+  (apply #'format stream "~d ~d ~d ~d" (vector->list (var-vector cv))))
38 37
 
39 38
 ;; вспомогательный макрос для консервативных переменных
40 39
 
... ...
@@ -156,11 +155,14 @@
156 156
 
157 157
 ;; операция minmod - особый случай
158 158
 
159
-(defmethod minmod ((a number) (b number))
160
-  (if (< (abs a) (abs b)) a b))
159
+(defmethod minmod ((type (eql 'number)) &rest params)
160
+  (cond ((andmap #'plusp params) (apply #'min params))
161
+	((andmap #'minusp params) (apply #'max params))
162
+	(t 0)))
161 163
 
162
-(defmethod minmod ((a conservative-variables) (b conservative-variables))
163
-  (conservative-variables-operator (list a b) #'minmod))
164
+(defmethod minmod ((type (eql 'conservative-variables)) &rest params)
165
+  (conservative-variables-operator params (curry #'minmod 'number)))
166
+
167
+(defmethod minmod ((type (eql 'primitive-variables)) &rest params)
168
+  (primitive-variables-operator params (curry #'minmod 'number)))
164 169
 
165
-(defmethod minmod ((a primitive-variables) (b primitive-variables))
166
-  (primitive-variables-operator (list a b) #'minmod))