Browse code

Добавил расчёт невязки и условие останова на по ней.

Dmitrii Kashin authored on 26/05/2015 00:41:14
Showing 9 changed files
... ...
@@ -76,14 +76,17 @@
76 76
 		 :number 'undef))
77 77
 
78 78
 (defmethod clone ((cell cell2d-cv))
79
-  (make-instance 'cell2d-cv
80
-		 :coord-vector (coord-vector cell)
81
-		 :size-vector (size-vector cell)
82
-		 :P (P cell)
83
-		 :phi (phi cell)
84
-		 :ug (ug cell)
85
-		 :ul (ul cell)
86
-		 ))
79
+  (let ((duplicate (make-instance 'cell2d-cv
80
+				  :coord-vector (coord-vector cell)
81
+				  :size-vector (size-vector cell)
82
+				  :P (P cell)
83
+				  :phi (phi cell)
84
+				  :ug (ug cell)
85
+				  :ul (ul cell)
86
+				  )))
87
+    (setf (convergence duplicate)
88
+	  (convergence cell))
89
+    duplicate))
87 90
 
88 91
 
89 92
 ;; строго говоря, не правильно добавлять градиенты в cell2d-cv. Надо
... ...
@@ -10,7 +10,7 @@
10 10
   :long-description "Данный пакет содержит реализацию решения задачи
11 11
   Фролова в примитивных переменных."
12 12
   :components ((:file "package-clf")
13
-	       (:file "framework")
13
+	       (:file "framework" :depends-on ("constants"))
14 14
 	       (:file "constants")
15 15
 	       (:file "variables-primitive" :depends-on ("constants" "framework"))
16 16
 	       (:file "cell-coordinates" :depends-on ("framework"))
... ...
@@ -15,7 +15,8 @@
15 15
 ;; конфигурация программы
16 16
 
17 17
 (defparameter *number-of-primitive-variables* 4)
18
-(defparameter *format-variables* (string-join (loop repeat *number-of-primitive-variables* collect "~d")))
18
+;;(defparameter *format-variables* (string-join (loop repeat *number-of-primitive-variables* collect "~d")))
19
+(defparameter *format-variables* "~d ~d ~d ~d")
19 20
 
20 21
 (defparameter *number-of-conservative-variables* 4)
21 22
 
... ...
@@ -310,3 +310,10 @@
310 310
 
311 311
 (defgeneric minmod (type &rest params)
312 312
   (:documentation "Вычисляет минимум по модулю"))
313
+
314
+(defgeneric count-convergence (old new)
315
+  (:documentation "Вычисляет невязку по старым и новым
316
+  параметрам. Параметры, вообще говоря, могут быть как ячейками, так и
317
+  сетками. Если параметры являются ячейками, то метод должен
318
+  возвращать число. Если сетками, то вычисляется невязка в каждой
319
+  ячейке, после чего она записывается в новую сетку."))
... ...
@@ -22,6 +22,10 @@
22 22
    ;; массив ячеек
23 23
    (cells :initform 'undef
24 24
 	  :accessor cells)
25
+   ;; максимальная невязка на сетке
26
+   (convergence :initform 0.0
27
+		:initarg :convergence
28
+		:accessor convergence)
25 29
    ;; граничные условия
26 30
    (left-boundary-condition :initarg :left-boundary-condition
27 31
 			    :initform (error "не задано левое ГУ")
... ...
@@ -275,6 +279,7 @@
275 275
   ;; информация о текущем состоянии
276 276
   (format t "# Iterations: ~a~%" (iteration grid))
277 277
   (format t "# Time: ~a~%" (grid-time grid))
278
+  (format t "# Convergence: ~a~%" (convergence grid))
278 279
   ;; заголовки
279 280
   (format t "NX NY ")
280 281
   (print-headers t (cell grid (index 0 0)))
... ...
@@ -302,7 +307,8 @@
302 302
 				:bottom-boundary-condition (bottom-boundary-condition grid)
303 303
 				:top-boundary-condition (top-boundary-condition grid)
304 304
 				:iteration (iteration grid)
305
-				:time (grid-time grid))))
305
+				:time (grid-time grid)
306
+				:convergence (convergence grid))))
306 307
     ;; скопируем все ячейки
307 308
     (loop for i from -1 to (1+ (cells-max-x grid)) do
308 309
 	 (loop for j from -1 to (1+ (cells-max-y grid)) do
... ...
@@ -1,6 +1,7 @@
1 1
 (in-package :clf)
2 2
 
3
-(defun march (grid &key output-every-step output-every-time finish-time finish-step
3
+(defun march (grid &key output-every-step output-every-time
4
+		     finish-time finish-step finish-convergence
4 5
 		     print-out-border suppress-export
5 6
 		     (method :kurganov-tadmor-explicit)
6 7
 		     (cu nil))
... ...
@@ -16,6 +17,8 @@
16 16
     (format t "Остановка после ~d итераций.~%" finish-step))
17 17
   (when finish-time
18 18
     (format t "Остановка после ~a секунд.~%" finish-time))
19
+  (when finish-convergence
20
+    (format t "Остановка при невязке ~a~%" finish-convergence))
19 21
 
20 22
   (let ((stop-flag nil)
21 23
 	(newgrid (clone grid)))
... ...
@@ -63,11 +66,24 @@
63 63
 	     (setf stop-flag t))
64 64
 	   
65 65
 	   ;; шаг
66
-	   (format t "Iter: ~d, Tau: ~d, Time: ~d~%" iterations tau time)
66
+	   (format t "Iter: ~d, Tau: ~d, Time: ~d" iterations tau time)
67
+	   (when finish-convergence (format t ", Convergence: ~d" (convergence grid)))
68
+	   (newline)
67 69
 	   (make-step grid newgrid tau method)
70
+	   
68 71
 	   ;; заботимся о выставлении параметров сетки
69 72
 	   (incf (iteration newgrid))
70 73
 	   (incf (grid-time newgrid) tau)
74
+
75
+	   ;; проверяем, достигли ли мы необходимой для остановки расчёта невязки
76
+	   (count-convergence grid newgrid)
77
+	   (when (and finish-convergence
78
+		      (< (convergence newgrid) finish-convergence))
79
+	     (format t "Достигнута сходимость по невзяке. Остановка.~%")
80
+	     (format t "Невязка сетки: ~a~%Требуемая невязка: ~a~%" (convergence newgrid) finish-convergence)
81
+	     (setf print-flag t)
82
+	     (setf stop-flag t))
83
+
71 84
 	   ;; печать
72 85
 	   (when (and (not suppress-export)
73 86
 		      print-flag)
... ...
@@ -14,18 +14,18 @@
14 14
     (send-output-to-file filename (print-grid grid :print-out-border print-out-border))
15 15
     (format t "OK~%")))
16 16
 
17
-(defun export-nodes (grid)
18
-  (let ((filename (make-nodes-filename grid)))
19
-    (format t "Export nodes into file ~s... " filename)
20
-    (send-output-to-file filename (print-nodes grid))
21
-    (format t "OK~%")))
22
-
23
-(defun export-grid&nodes (grid)
24
-  (export-grid grid :print-out-border t)
25
-  (export-nodes grid))
17
+;;(defun export-nodes (grid)
18
+;;  (let ((filename (make-nodes-filename grid)))
19
+;;    (format t "Export nodes into file ~s... " filename)
20
+;;    (send-output-to-file filename (print-nodes grid))
21
+;;    (format t "OK~%")))
22
+;;
23
+;;(defun export-grid&nodes (grid)
24
+;;  (export-grid grid :print-out-border t)
25
+;;  (export-nodes grid))
26 26
 
27 27
 (defun gridbg (text grid)
28 28
   (when *debug-grids-in-make-step*
29
-    (format t "~a: ~%" text)
29
+    (format t "ttt ~a: ~%" text)
30 30
     (print-grid grid :print-out-border nil)
31 31
     (format t "----~%~%")))
... ...
@@ -461,6 +461,25 @@
461 461
     ;;newgrid ; в будущем планирую возвращать сетку, а не писваивать значения. Может зря.
462 462
     ))
463 463
 
464
+;; вычисление невязки
465
+(defmethod count-convergence ((oldcell cell2d-cv) (newcell cell2d-cv))
466
+  (let ((oldcell (conservative->primitive oldcell))
467
+	(newcell (conservative->primitive newcell)))
468
+    (apply #'max (vector->list (var-vector (pvar- newcell oldcell))))))
469
+
470
+(defmethod count-convergence ((oldgrid grid2d-cv) (newgrid grid2d-cv))
471
+  (setf (convergence newgrid) 0.0d0)
472
+  (loop for j from (min-y-index oldgrid) to (max-y-index oldgrid) do
473
+       (loop for i from (min-x-index oldgrid) to (max-x-index oldgrid)
474
+	  for index = (index i j)
475
+	  for oldcell = (cell oldgrid index)
476
+	  for newcell = (cell newgrid index)
477
+	  for conv = (count-convergence oldcell newcell)
478
+	  do
479
+	    (setf (convergence newcell) conv)
480
+	    (when (> conv (convergence newgrid))
481
+	      (setf (convergence newgrid) conv)))))
482
+
464 483
 ;; солвер для неявной схемы курганова-тадмора
465 484
 ;(defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number) (method (eql :kurganov-tadmor)))
466 485
 ;  (labels
... ...
@@ -85,7 +85,64 @@
85 85
       (gridbg "4 (boundaries)" grid)
86 86
       (gridbg "start-grid" grid)
87 87
       (let ((newgrid (march grid
88
-			    :finish-time 0.05
88
+			    :finish-time 0.005
89
+			    :output-every-time 0.0001
90
+			    :cu cu
91
+			    )))
92
+	(gridbg "final-grid" newgrid)
93
+	newgrid
94
+	)
95
+      )))
96
+
97
+(defun test-kt-split-1/convergence (&key (cells-nx 100) (cells-ny 1) (cu 0.001)
98
+				      (output-dir #P"~/prog/clfpv/test-kt-split-1-convergence/"))
99
+  "Тест на стоячем скачке уплотнения. В данном тесте значение phi
100
+близко к 1, а это значит, что решение должно быть почти таким же, как
101
+в газе. Отличие от test-kt-split-1 заключается в том, что здесь я
102
+тестирую остановку по невязке."
103
+  (ensure-directories-exist output-dir)
104
+  (let ((*DEFAULT-PATHNAME-DEFAULTS* output-dir))
105
+    (let ((grid (make-instance 'grid2d-cv
106
+			       :cells-amount-x cells-nx
107
+			       :cells-amount-y cells-ny
108
+			       :left-boundary-condition 'zeroe
109
+			       :right-boundary-condition 'zeroe
110
+			       :bottom-boundary-condition 'zeroe
111
+			       :top-boundary-condition 'zeroe)))
112
+      (gridbg "1 (instance)" grid)
113
+      (init-grid grid
114
+		 :left-border-coord 0.0d0
115
+		 :right-border-coord 1.0d0
116
+		 :top-border-line (constant-function 1.0d0)
117
+		 :bottom-border-line (constant-function 0.0d0)
118
+		 :thickening-coef-x 1
119
+		 :thickening-coef-y 1)
120
+      (gridbg "2 (init-grid)" grid)
121
+      (let* ((left-vars (make-instance 'primitive-variables
122
+				       :phi 0.999999
123
+				       :rg 1.0
124
+				       :ug 400.0
125
+				       :ul 400.0))
126
+	     (right-vars (make-instance 'primitive-variables
127
+					:phi (phi left-vars)
128
+					:ug (/ (* *sos* *sos*)
129
+					       (* *kappa* (ug left-vars)))
130
+					:ul (/ (* *sos* *sos*)
131
+					       (* *kappa* (ul left-vars)))
132
+					:rg (/ (* (rg left-vars) (ug left-vars))
133
+					       ;; на ug справа
134
+					       (/ (* *sos* *sos*)
135
+						  (* *kappa* (ug left-vars)))))))
136
+	(set-splitted-field-2d-x grid
137
+				 :left-vars left-vars
138
+				 :right-vars right-vars))
139
+      (gridbg "3 (set-field)" grid)
140
+      (update-boundary-cells grid)
141
+      (gridbg "4 (boundaries)" grid)
142
+      (gridbg "start-grid" grid)
143
+      (let ((newgrid (march grid
144
+			    :finish-time 0.005
145
+			    :finish-convergence 1.0d-6
89 146
 			    :output-every-time 0.0001
90 147
 			    :cu cu
91 148
 			    )))