Browse code

многочисленные исправления ошибок

(запустил тест и прошёлся по ругани интерпретатора)

Dmitrii Kashin authored on 11/05/2015 05:43:57
Showing 9 changed files
... ...
@@ -23,16 +23,12 @@
23 23
 	       (:file "grid2d-conservative" :depends-on ("cells-primitive" "cell-index"))
24 24
 	       (:file "field2d-conservative" :depends-on ("grid2d-conservative"))
25 25
 	       (:file "boundary-conditions" :depends-on ("grid2d-primitive"))
26
-	       
27
-;;	       (:file "solver" :depends-on ("package" "field"))
28
-;;	       (:file "solver2d" :depends-on ("package" "field2d"))
29
-;;	       (:file "march2d" :depends-on ("package" "solver" "grid" "field"))
30
-
31 26
 	       (:file "solver2d-kurganov-tadmor" :depends-on ("grid2d-conservative"))
32
-	       (:file "test2d-conservative" :depends-on ("field2d-conservative" "boundary-conditions"))
27
+	       (:file "march" :depends-on ("grid2d-conservative" "output-functions"))
28
+	       (:file "test2d-conservative" :depends-on ("field2d-conservative" "boundary-conditions" "march"))
33 29
 	       (:file "output-functions")
34 30
 	       
35
-	       (:file "march" :depends-on ("grid2d-conservative" "output-functions"))
31
+	       
36 32
 	       ))
37 33
   
38 34
 
... ...
@@ -1,50 +1,35 @@
1 1
 (in-package :clf)
2 2
 
3
-;; монотонное поле
4
-
5
-(defmacro set-whole-field-2d (grid pv)
6
-  `(loop for i from 0 to (cells-max-x grid) do
7
-	(loop for j from 0 to (cells-max-y grid) do
8
-	     (setf (var-vector (cell grid (index i j)))
9
-		   (var-vector (primitive->conservative ,pv))))))
10
-
11
-(defun set-constant-field-2d (grid)
3
+(defun set-constant-field-2d (grid pv)
12 4
   "Задаёт в сетке `grid' постоянное поле"
13
-  )
14
-
15
-
16
-;; этот макрос принимает левые и правые консервативные переменные, а
17
-;; потому задаёт поле, линейно изменяющееся от левых параметров к
18
-;; правым.
19
-
20
-;; todo: надо бы сделать проверки на случай не задания left-vars и right-vars
21
-
22
-(defmacro set-monotonous-field-2d-x (grid &key left-vars right-vars)
23
-  "Задаёт в сетке `grid' линейно изменяющееся поле"
24
-  `(let ((pv-grad (pvar/ (pvar- ,right-vars ,left-vars)
5
+  (loop for i from 0 to (cells-max-x grid) do
6
+       (loop for j from 0 to (cells-max-y grid) do
7
+	    (setf (cell grid (index i j))
8
+		  (primitive->conservative pv)))))
9
+
10
+(defun set-monotonous-field-2d-x (grid &key left-vars right-vars)
11
+  "Задаёт в сетке `grid' поле, линейно изменяющееся от `left-vars' до `right-vars'"
12
+  (let ((pv-grad (pvar/ (pvar- right-vars left-vars)
25 13
 			 (cells-amount-x grid))))
26 14
      (loop for i from 0 to (cells-max-x grid) do
27
-	  (let ((pv (pvar+ ,left-vars
15
+	  (let ((pv (pvar+ left-vars
28 16
 			   (pvar* (+ i 0.5) pv-grad))))
29 17
 	    (loop for j from 0 to (cells-max-y grid) do
30 18
 	       (setf (cell grid (index i j))
31 19
 		     (primitive->conservative pv)))))))
32 20
 
33
-
34
-;; этот макрос нужен для тестов на скачах
35
-
36
-(defmacro set-splitted-field-2d-x (grid &key left-vars right-vars)
37
-  `(progn
38
-     (format t "Задаём поле, разделённое по оси x~%")
39
-     (format t "Колонки: ~T") (print-headers t ,left-vars) (newline)
40
-     (format t "left:~T") (print-content t ,left-vars :with-headers nil) (newline)
41
-     (format t "right:~T") (print-content t ,right-vars :with-headers nil) (newline)
42
-     (let* ((mid-x-index (round (/ (cells-amount-x ,grid) 2)))
43
-	    (left-vars (primitive->conservative ,left-vars))
44
-	    (right-vars (primitive->conservative ,right-vars)))
45
-       (loop for i from 0 to (cells-max-x grid) do
46
-	    (loop for j from 0 to (cells-max-y grid) do
47
-		 (setf (cell grid (index i j))
48
-		       (if (< i mid-x-index)
49
-			   left-vars
50
-			   right-vars)))))))
21
+(defun set-splitted-field-2d-x (grid &key left-vars right-vars)
22
+  "Задаёт в сетке `grid' поле, до половины заполненное `left-vars', а после - `right-vars'"
23
+  (format t "Задаём поле, разделённое по оси x~%")
24
+  (format t "Колонки: ~T") (print-headers t left-vars) (newline)
25
+  (format t "left:~T") (print-content t left-vars :with-headers nil) (newline)
26
+  (format t "right:~T") (print-content t right-vars :with-headers nil) (newline)
27
+  (let* ((mid-x-index (round (/ (cells-amount-x grid) 2)))
28
+	 (left-vars (primitive->conservative left-vars))
29
+	 (right-vars (primitive->conservative right-vars)))
30
+    (loop for i from 0 to (cells-max-x grid) do
31
+	 (loop for j from 0 to (cells-max-y grid) do
32
+	      (setf (cell grid (index i j))
33
+		    (if (< i mid-x-index)
34
+			left-vars
35
+			right-vars))))))
... ...
@@ -283,7 +283,13 @@
283 283
   (:documentation "Возвращает вектор источниковых членов, рассчитанный
284 284
   по переменным `vars'"))
285 285
 
286
-;; --------------------------------------------------
287
-
288 286
 (defgeneric update-boundary-cells (grid)
289
-  (:documentation "Set boundary cells according to boundary conditions"))
287
+  (:documentation "Устанавливает в сетке `grid' заграничные ячейки в
288
+  соответствии с её граничными условиями"))
289
+
290
+(defgeneric minmod (a b)
291
+  (:documentation "Вычисляет минимум по модулю"))
292
+
293
+
294
+
295
+
... ...
@@ -55,6 +55,9 @@
55 55
 ;; note: у меня такое чувство, будто мутаторами я пользоваться не буду
56 56
 
57 57
 (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)))
58 61
   (aref (cells grid) (1+ (x-index index)) (1+ (y-index index))))
59 62
 
60 63
 (defmethod (setf cell) ((new cell2d-cv) (grid grid2d-cv) (index cell-index))
... ...
@@ -13,6 +13,10 @@
13 13
 (push "/usr/share/common-lisp/source/alexandria/" asdf:*central-registry*)
14 14
 ;(push "/home/freehck/quicklisp/dists/quicklisp/software/gsll-20140211-git" asdf:*central-registry*)
15 15
 
16
+; если у меня руки дойдёт вычислять производные автоматически
17
+;(push "/home/freehck/prog/maxima-5.35.1/src/" asdf:*central-registry*)
18
+;(asdf:load-system :maxima)
19
+
16 20
 (asdf::load-asd "clf.asd")
17 21
 ;(asdf::load-asd "alexandria.asd")
18 22
 
... ...
@@ -17,8 +17,6 @@
17 17
     (format nil "Остановка после ~d итераций.~%" finish-step))
18 18
   (when finish-time
19 19
     (format t "Остановка после ~a секунд.~%" finish-time))
20
-  (when (not cu)
21
-    (error "Число CFL не задано"))
22 20
 
23 21
   (let ((stop-flag nil)
24 22
 	(newgrid (clone grid)))
... ...
@@ -174,27 +174,36 @@
174 174
 		   :var-vector (vector 0 sg 0 (- sg)))))
175 175
 
176 176
 ;; градиенты
177
-(defun minmod (a b)
178
-  (if (< (abs a) (abs b)) a b))
179
-
180 177
 (defun grad-x (grid indx order)
181 178
   (let* ((i (x-index indx))
182 179
 	 (j (y-index indx))
180
+	 (outborder-left (< i 0))
181
+	 (outborder-right (>= i (cells-amount-x grid)))
183 182
 	 (mid-cell (cell grid indx))
184
-	 (left-cell (cell grid (index (1- i) j)))
185
-	 (right-cell (cell grid (index (1+ i) j)))
186
-	 (DL (- (x-coord right-cell) (x-coord mid-cell)))
187
-	 (DR (- (x-coord mid-cell) (x-coord left-cell)))
188
-	 (DUL (- right-cell mid-cell))
189
-	 (DUR (- mid-cell left-cell)))
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))))
190 195
     (case order
191
-      ((1) (minmod (/ DUL DL)
192
-		   (/ DUR DR)))
193
-      ((2) (/ (+ (/ (* DUL DR)
194
-		    DL)
195
-		 (/ (* DUR DL)
196
-		    DR))
197
-	      (+ DL DR)))
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)))
198 207
       (otherwise (error "`grad' умеет считать производные лишь до 2го порядка включительно")))))
199 208
 
200 209
 (defun grad1-x (grid indx)
... ...
@@ -231,14 +240,15 @@
231 231
 					       u3)
232 232
 					    (* 2 ul u4)))))))
233 233
 
234
+;; todo: а правильно ли я стреляю? линейны-то примитивные переменные
234 235
 ;; выстрелы на границу
235 236
 (defun shot (grid indx direction)
236 237
   (let* ((cell (cell grid indx))
237 238
 	 (hx (x-size cell))
238
-	 (gx (grad1-x cell indx)))
239
+	 (gx (grad1-x grid indx)))
239 240
     (case direction
240
-      ((left) (- cell (* 0.5 hx gx)))
241
-      ((right) (+ cell (* 0.5 hx gx)))
241
+      ((left) (cvar- cell (cvar* 0.5 hx gx)))
242
+      ((right) (cvar+ cell (cvar* 0.5 hx gx)))
242 243
       (otherwise (error 'shot "Неизвестное направление для стрельбы")))))
243 244
 
244 245
 (defun left-shot (grid indx)
... ...
@@ -257,9 +267,13 @@
257 257
 	 (j (y-index indx))
258 258
 	 (l-indx (index (1- i) j))
259 259
 	 (r-indx (index (1+ i) j))
260
+	 (outborder-left (< i 0))
261
+	 (outborder-right (>= i (cells-amount-x grid)))
260 262
 	 (mid-cell (cell grid indx))
261
-	 (left-cell (cell grid l-indx))
262
-	 (right-cell (cell grid r-indx)))
263
+	 (left-cell (unless outborder-left
264
+		      (cell grid l-indx)))
265
+	 (right-cell (unless outborder-right
266
+		       (cell grid r-indx))))
263 267
     (case direction
264 268
       ((right) (max (radius (right-shot grid indx) 'A)
265 269
 		    (radius (left-shot grid r-indx) 'A)
... ...
@@ -310,7 +324,6 @@
310 310
     (border-flux grid indx tau :left)))
311 311
 
312 312
 ;; солвер
313
-
314 313
 (defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number))
315 314
   (labels
316 315
       ((predictor (grid tau)
... ...
@@ -324,24 +337,45 @@
324 324
 			    (u (cell grid (index i j)))
325 325
 			    (hx (x-size u))
326 326
 			    (s (source u)))
327
-		       (flet ((eq-system-non-linear (x1 x2 x3 x4)
328
-				(let* ((x (make-instance 'conservative-variables
329
-							 :var-vector (vector x1 x2 x3 x4)))
330
-				       (result (cvar- (cvar+ (cvar/ (cvar- x u)
331
-								    tau)
332
-							     (/ dF hx))
333
-						      (* 0.5 (cvar+ s (source x))))))
334
-				  (let-conservative result
335
-				    (values u1 u2 u3 u4)))))
327
+		       (labels (;; функция, для которой мы будем искать нули методом Ньютона
328
+				(eq-system (x1 x2 x3 x4)
329
+				  (let* ((x (make-instance 'conservative-variables
330
+							   :var-vector (vector x1 x2 x3 x4)))
331
+					 (result (cvar- (cvar+ (cvar/ (cvar- x u)
332
+								      tau)
333
+							       (cvar/ dF hx))
334
+							(cvar* 0.5 (cvar+ s (source x))))))
335
+				    (let-conservative result
336
+				      (values u1 u2 u3 u4))))
337
+				;; частные производеные этой функции (якобиан)
338
+				(eq-system-df (x1 x2 x3 x4)
339
+				  (let* ((D (+ (/ x2 x1)
340
+					       (/ x4 x3)))
341
+					 (F (* 0.5 *rl* D (signum D)))
342
+					 (F1 (/ (* F x2) (* x1 x1)))
343
+					 (F2 (/ F x1))
344
+					 (F3 (/ (* F x4) (* x3 x3)))
345
+					 (F4 (/ F x3))
346
+					 (TT (/ 1 tau)))
347
+				      (values TT 0 0 0 ;;
348
+					      F1 (- TT F2) (- F3) F4 ;;
349
+					      0 0 TT 0 ;;
350
+					      (- F1) F2 F3 (- TT F4))))
351
+				;; собирательная функция (необходима для GSLL)
352
+				(eq-system-fdf (x1 x2 x3 x4)
353
+				  (multiple-value-bind (v1 v2 v3 v4)
354
+				      (eq-system x1 x2 x3 x4)
355
+				    (multiple-value-bind (j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)
356
+					(eq-system-df x1 x2 x3 x4)
357
+				      (values v1 v2 v3 v4 j1 j2 j3 j4 j5 j6 j7 j8 j9 j10 j11 j12 j13 j14 j15 j16)))))
336 358
 			 (let-conservative u
337 359
 			   (let ((max-iter *multi-root-iter*)
338 360
 				 (solver (make-multi-dimensional-root-solver-f 
339
-					  method #'eq-system-non-linear
361
+					  +newton-mfdfsolver+
362
+					  (list #'eq-system #'eq-system-df #'eq-system-fdf)
340 363
 					  (grid:make-foreign-array
341 364
 					   'double-float :initial-contents (list u1 u2 u3 u4)))))
342 365
 			     (loop for iter from 0
343
-				;;for fnval = nil then (function-value solver)
344
-				;;for argval = nil then (solution solver)
345 366
 				while (if (< iter max-iter)
346 367
 					  (or (zerop iter)
347 368
 					      (not (multiroot-test-residual solver *multi-root-est-error*))))
... ...
@@ -364,14 +398,15 @@
364 364
 		   for dF+ = (border-flux midgrid index tau :right)
365 365
 		   for dF- = (border-flux midgrid index tau :left) then dF+
366 366
 		   do
367
-		     (let ((S (cvar* 0.5 (cvar+ (source grid index)
368
-						(source predgrid index))))
369
-			   (u (cell grid (index i j)))
370
-			   (hx (x-size u))
371
-			   (dF (cvar- dF+ dF-)))
367
+		     (let* ((u (cell grid index))
368
+			    (S (cvar* 0.5 (cvar+ (source u)
369
+						 (source (cell predgrid index)))))
370
+			    (hx (x-size u))
371
+			    (dF (cvar- dF+ dF-)))
372 372
 		       (setf (cell newgrid index)
373 373
 			     (cvar- (cvar+ u (cvar* tau S))
374 374
 				    (cvar/ (cvar* tau dF)
375 375
 					   hx))))))
376 376
 	   newgrid)))
377 377
     (corrector grid (predictor grid tau) tau)))
378
+
... ...
@@ -18,13 +18,12 @@
18 18
 		 :bottom-border-line (constant-function 0)
19 19
 		 :thickening-coef-x 1
20 20
 		 :thickening-coef-y 1)
21
-      (set-whole-field-2d grid
22
-			  (make-instance 'conservative-variables
23
-					 :phi 0.999999
24
-					 :P 100000.0
25
-					 :ug 25.0
26
-					 :ul 25.0))
27
-      (set-constant-field-2d grid)
21
+      (set-constant-field-2d grid
22
+			     (make-instance 'conservative-variables
23
+					    :phi 0.999999
24
+					    :P 100000.0
25
+					    :ug 25.0
26
+					    :ul 25.0))
28 27
       (update-boundary-cells grid)
29 28
       (march grid :finish-step -1)
30 29
       grid
... ...
@@ -109,3 +108,38 @@
109 109
       (march grid :finish-step -1)
110 110
       grid
111 111
       )))
112
+
113
+
114
+(defun test-kt-constant (&key (cells-nx 20) (cells-ny 1) (output-dir #P"~/prog/clfpv/test-KT-constant/"))
115
+  "Тест на постоянном поле"
116
+  (ensure-directories-exist output-dir)
117
+  (let ((*DEFAULT-PATHNAME-DEFAULTS* output-dir))
118
+    (let ((grid (make-instance 'grid2d-cv
119
+			       :cells-amount-x cells-nx
120
+			       :cells-amount-y cells-ny
121
+			       :left-boundary-condition 'zeroe
122
+			       :right-boundary-condition 'zeroe
123
+			       :bottom-boundary-condition 'zeroe
124
+			       :top-boundary-condition 'zeroe)))
125
+      (init-grid grid
126
+		 :left-border-coord 0
127
+		 :right-border-coord 1
128
+		 :top-border-line (constant-function 1)
129
+		 :bottom-border-line (constant-function 0)
130
+		 :thickening-coef-x 1
131
+		 :thickening-coef-y 1)
132
+      (set-constant-field-2d grid
133
+			     (make-instance 'conservative-variables
134
+					    :phi 0.999999
135
+					    :P 100000.0
136
+					    :ug 25.0
137
+					    :ul 25.0))
138
+      (update-boundary-cells grid)
139
+      (march grid :finish-step 1)
140
+      grid
141
+      )))
142
+
143
+
144
+
145
+
146
+    
... ...
@@ -114,8 +114,8 @@
114 114
 		     :rg rg
115 115
 		     :ug ug
116 116
 		     :ul ul))))
117
-;; операции, которые можно производить над консервативными переменными
118 117
 
118
+;; операции, которые можно производить над консервативными переменными
119 119
 (defun number->conservative (num)
120 120
   (make-instance 'conservative-variables :number num))
121 121
 
... ...
@@ -129,7 +129,7 @@
129 129
     `(progn
130 130
        (when (null ,cv-list) (error "conservative-variables-operator нельзя применять к пустому списку"))
131 131
        (let ((,cv-list-sym (mapcar #'any->conservative ,cv-list)))
132
-	 (make-instance 'primitive-variables
132
+	 (make-instance 'conservative-variables
133 133
 			:var-vector (apply #'map 'vector ,func
134 134
 					   (mapcar #'var-vector ,cv-list-sym)))))))
135 135
 
... ...
@@ -153,3 +153,11 @@
153 153
 
154 154
 (defun cvar-signum (cv)
155 155
   (conservative-variables-operator (list cv) #'signum))
156
+
157
+;; операция minmod - особый случай
158
+
159
+(defmethod minmod ((a number) (b number))
160
+  (if (< (abs a) (abs b)) a b))
161
+
162
+(defmethod minmod ((a conservative-variables) (b conservative-variables))
163
+  (conservative-variables-operator (list a b) #'minmod))