Browse code

Доделываю Курганова-Тадмора. Почти все основные функции + сделал диспетчеризатор для поиска собственных чисел.

Dmitrii Kashin authored on 07/05/2015 12:29:45
Showing 2 changed files
... ...
@@ -6,3 +6,4 @@
6 6
 	:gsll)
7 7
   ;; следующие символы конфликтовали между alexandria и gsll
8 8
   (:shadow mean variance median standard-deviation factorial))
9
+
... ...
@@ -1,10 +1,5 @@
1 1
 (in-package :clf)
2 2
 
3
-;; определение шага по времени
4
-
5
-(defmethod max-timestep-constant ((grid grid2d-cv) (method (eql :lax-frederichs)) cu)
6
-  5.e-8)
7
-
8 3
 ;; функция-добавок к давлению
9 4
 (defun sigma (phi)
10 5
   phi *sigma*)
... ...
@@ -13,131 +8,108 @@
13 13
 (defun Dsigma (phi)
14 14
   phi 0)
15 15
 
16
-
17
-(defun eigen-values-by-cell (cell)
16
+;; собственные числа ячейки
17
+;; возвращает одно из 3 значений, в зависимости от выбранного type: 'A, 'S или 'Afull
18
+(defun eigen-values-by-cell (cell type)
18 19
   (let* ((phi (phi cell))
19 20
 	 (rg (rg cell))
20 21
 	 (ug (ug cell))
21 22
 	 (ul (ul cell))
22 23
 	 (rl *rl*)
23
-	 (cg (/ *Cg* (sqrt *kappa*)))
24
-	 (cl (sqrt (/ (+ (* rg cg cg)
25
-			 (* phi (- (sigma phi)
26
-				   rl ul ul)))
27
-		      (* phi rl))))
28
-	 (eigen-values/A ; !!!
29
-	  (list (abs (- ug cg))
30
-		(abs (+ ug cg))
31
-		(abs (+ ul cl))
32
-		(abs (- ul cl))))
33
-	 (eigen-values/sources ; !!!
34
-	  (list (* -1
35
-		   (abs (- ug ul))
36
-		   (+ (/ *rl*
37
-			 (* phi rg))
38
-		      (/ 1
39
-			 (- 1 phi))))))
40
-	 (suggestion (* 5 cg))
41
-	 (left-init (- (min ul (- ug cg) (+ ug cg))
42
-		       suggestion))
43
-	 (right-init (+ (max ul (- ug cg) (+ ug cg))
44
-			suggestion))
45
-	 
46
-	 (X (- (sigma phi)
47
-	       (* (1- phi)
48
-		  (Dsigma phi))))
49
-	 (Z (* (1- phi) rg))
50
-	 (R (+ (* Z cg cg)
51
-		(* X phi)))
52
-	 (c1 (* phi rl))
53
-	 (c2 (* -2 c1 (+ ug ul)))
54
-	 (c3 (- (* c1 (+ (expt (+ ug ul) 2)
55
-			 (* 2 ug ul)
56
-			 (* -1 cg cg)))
57
-		R))
58
-	 (c4 (- (* 2 ug R)
59
-		(* 2 c1 ul
60
-		   (- (* (+ ug ul) ug)
61
-		      (* cg cg)))))
62
-	 (c5 (- (* X cg cg phi)
63
-		(* R ug ug)
64
-		(* c1 ul ul (- (* cg cg)
65
-			       (* ug ug))))))
66
-    (flet ((frolov-equation (x)
67
-	     (- (* phi
68
-		   (- (* rl (expt (- x ul) 2))
69
-		      (- (sigma phi)
70
-			 (* (- 1 phi) (Dsigma phi))))
71
-		   (- (expt (- x ug) 2)
72
-		      (expt cg 2)))
73
-		(* (expt cg 2)
74
-		   (- 1 phi)
75
-		   rg
76
-		   (expt (- x ug) 2))))
77
-	   (frolov-equation-derivative (x)
78
-	     (- (+ (* 2 phi rl
79
-		      (- x ul)
80
-		      (- (expt (- x ug) 2)
81
-			 (expt cg 2)))
82
-		   (* 2 phi
83
-		      (- x ug)
84
-		      (- (* rl (expt (- x ul) 2))
85
-			 (- (sigma phi)
86
-			    (* (- 1 phi) (Dsigma phi))))))
87
-		(* 2 rg
88
-		   (expt cg 2)
89
-		   (- 1 phi)
90
-		   (- x ug))))
91
-	   (frolov-equation-fdf (x)
92
-	     (values (frolov-equation x)
93
-		     (frolov-equation-derivative x)))
94
-	   (find-root (&optional initial (method +newton-fdfsolver+))
95
-	     (let* ((max-iter *newton-max-iter*)
96
-		    (solver (make-one-dimensional-root-solver-fdf
97
-			     method 'frolov-equation 'frolov-equation-derivative 'frolov-equation-fdf initial)))
98
-	       (loop for iter from 0
99
-		  for oldroot = initial then root
100
-		  for root = (progn (iterate solver) (solution solver))
101
-		  while (if (< iter max-iter)
102
-			    (not (root-test-delta root oldroot *newton-est-error* 0.0d0))
103
-			    (error "Метод Ньютона превысил допустимое число итераций."))
104
-		  ;;для отдалки
105
-;		  do (format t "~d~6t~g ~34t~10,6g~&"
106
-;			     iter root (- root oldroot))
107
-		  finally (progn
108
-;			    (format t "~d~6t~g ~34t~10,6g~&~%"
109
-;				    iter root (- root oldroot))
110
-			    (return root))))))
111
-
112
-    (let ((eigen-values/Afull ; !!!
113
-	   (progn
114
-	     ;;(format t "rg: ~d~%" rg)
115
-	     ;;(format t "phi: ~d~%" phi)
116
-	     ;;(format t "ul: ~d~%" ul)
117
-	     ;;(format t "ug: ~d~%" ug)
118
-	     ;;(format t "rl: ~d~%" rl)
119
-	     ;;(format t "cg: ~d~%--~%~%" cg)
120
-	     ;;(format t "X: ~d~%" X)
121
-	     ;;(format t "Z: ~d~%" Z)
122
-	     ;;(format t "R: ~d~%" R)
123
-	     (let* ((root-1 (find-root left-init))
124
-		    (root-2 (find-root right-init))
125
-		    (coefficients (list c1 c2 c3 c4 c5))
126
-		    (quadro-coeffs (divide-polynomial coefficients root-1 root-2)))
127
-	       (multiple-value-bind (root-3 root-4)
128
-		   (apply #'solve-quadratic-complex
129
-			  (mapcar (curryr #'coerce 'double-float) quadro-coeffs))
130
-		 ;;(format t "root: ~d~%" root-1)
131
-		 ;;(format t "root: ~d~%" root-2)
132
-		 ;;(format t "coefficients: ~a~%" coefficients)
133
-		 ;;(format t "quadro-coeffs: ~a~%" quadro-coeffs)
134
-		 ;;(format t "root: ~d~%" root-3)
135
-		 ;;(format t "root: ~d~%" root-4)
136
-		 (list root-1 root-2 root-3 root-4))))))
137
-      
138
-      (values eigen-values/A
139
-	      eigen-values/sources
140
-	      eigen-values/Afull)))))
24
+	 (cg (/ *Cg* (sqrt *kappa*))))
25
+      (case type
26
+
27
+	((A flux) ;; собственные числа якобиана потоковых членов
28
+	 (let ((cl (sqrt (/ (+ (* rg cg cg)
29
+			       (* phi (- (sigma phi)
30
+					 rl ul ul)))
31
+			    (* phi rl)))))
32
+	   (list (abs (- ug cg))
33
+		 (abs (+ ug cg))
34
+		 (abs (+ ul cl))
35
+		 (abs (- ul cl)))))
36
+
37
+	((S source) ;; собственные числа (одно) якобиана источниковых членов
38
+	 (list (- (* (abs (- ug ul))
39
+		     (+ (/ rl
40
+			   (* phi rg))
41
+			(/ 1
42
+			   (- 1 phi)))))))
43
+
44
+	((Afull) ;; собственные числа якобиана потоковых членов, объединённых с
45
+		 ;; источниками, обусловленными сужением трубки тока
46
+	 (labels ((frolov-equation (x)
47
+		    (- (* phi
48
+			  (- (* rl (expt (- x ul) 2))
49
+			     (- (sigma phi)
50
+				(* (- 1 phi) (Dsigma phi))))
51
+			  (- (expt (- x ug) 2)
52
+			     (expt cg 2)))
53
+		       (* (expt cg 2)
54
+			  (- 1 phi)
55
+			  rg
56
+			  (expt (- x ug) 2))))
57
+		  (frolov-equation-derivative (x)
58
+		    (- (+ (* 2 phi rl
59
+			     (- x ul)
60
+			     (- (expt (- x ug) 2)
61
+				(expt cg 2)))
62
+			  (* 2 phi
63
+			     (- x ug)
64
+			     (- (* rl (expt (- x ul) 2))
65
+				(- (sigma phi)
66
+				   (* (- 1 phi) (Dsigma phi))))))
67
+		       (* 2 rg
68
+			  (expt cg 2)
69
+			  (- 1 phi)
70
+			  (- x ug))))
71
+		  (frolov-equation-fdf (x)
72
+		    (values (frolov-equation x)
73
+			    (frolov-equation-derivative x)))
74
+		  (find-root (&optional initial (method +newton-fdfsolver+))
75
+		    (let* ((max-iter *newton-max-iter*)
76
+			   (solver (make-one-dimensional-root-solver-fdf
77
+				    method #'frolov-equation #'frolov-equation-derivative #'frolov-equation-fdf initial)))
78
+		      (loop for iter from 0
79
+			 for oldroot = initial then root
80
+			 for root = (progn (iterate solver) (solution solver))
81
+			 while (if (< iter max-iter)
82
+				   (not (root-test-delta root oldroot *newton-est-error* 0.0d0))
83
+				   (error "Метод Ньютона превысил допустимое число итераций."))
84
+			 finally (return root)))))
85
+	   (let* ((suggestion (* 5 cg))
86
+		  (left-init (- (min ul (- ug cg) (+ ug cg))
87
+				suggestion))
88
+		  (right-init (+ (max ul (- ug cg) (+ ug cg))
89
+				 suggestion))
90
+		  (root-1 (find-root left-init))
91
+		  (root-2 (find-root right-init))
92
+		  (X (- (sigma phi)
93
+			(* (1- phi)
94
+			   (Dsigma phi))))
95
+		  (Z (* (1- phi) rg))
96
+		  (R (+ (* Z cg cg)
97
+			(* X phi)))
98
+		  (c1 (* phi rl))
99
+		  (c2 (* -2 c1 (+ ug ul)))
100
+		  (c3 (- (* c1 (+ (expt (+ ug ul) 2)
101
+				  (* 2 ug ul)
102
+				  (* -1 cg cg)))
103
+			 R))
104
+		  (c4 (- (* 2 ug R)
105
+			 (* 2 c1 ul
106
+			    (- (* (+ ug ul) ug)
107
+			       (* cg cg)))))
108
+		  (c5 (- (* X cg cg phi)
109
+			 (* R ug ug)
110
+			 (* c1 ul ul (- (* cg cg)
111
+					(* ug ug)))))
112
+		  (coefficients (list c1 c2 c3 c4 c5))
113
+		  (quadro-coeffs (divide-polynomial coefficients root-1 root-2)))
114
+	     (multiple-value-bind (root-3 root-4)
115
+		 (apply #'solve-quadratic-complex
116
+			(mapcar (curryr #'coerce 'double-float) quadro-coeffs))
117
+	       (list root-1 root-2 root-3 root-4))))))))
141 118
 
142 119
 
143 120
 ;; чертовски важная штука, без неё reader macro #m не заработает нормально
... ...
@@ -166,15 +138,149 @@
166 166
 
167 167
 ;; небольшие тесты
168 168
 
169
-;(defparameter *eigen-cell* (make-instance 'conservative-variables
170
-;					  :phi 0.01
171
-;					  :ug 25
172
-;					  :ul 400
173
-;					  :rg 0.001
174
-;					  ))
175
-;
176
-;(eigen-values-by-cell *eigen-cell*)
177
-						
169
+;;(defparameter *eigen-cell* (make-instance 'conservative-variables
170
+;;					  :phi 0.01
171
+;;					  :ug 25
172
+;;					  :ul 400
173
+;;					  :rg 0.001
174
+;;					  ))
175
+;;
176
+;;(eigen-values-by-cell *eigen-cell* 'A)
177
+;;(eigen-values-by-cell *eigen-cell* 'S)
178
+;;(eigen-values-by-cell *eigen-cell* 'Afull)
179
+
180
+;; определение шага по времени
181
+
182
+(defmethod max-timestep-constant ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
183
+  5.e-8)
184
+
185
+;; к сожалени, пока нет описания, как это делать
178 186
 (defmethod max-timestep/cu ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
179 187
   (noop))
180 188
 
189
+;; потоковый член в ячейке.
190
+(defmethod flux ((vars conservative-variables))
191
+  (let-conservative vars
192
+    (make-instance 'conservative-variables
193
+		   :var-vector (vector u2
194
+				       (+ (/ (* u2 u2)
195
+					     u1)
196
+					  (/ (* *Cg* *Cg* u1)
197
+					     *kappa*))
198
+				       u4
199
+				       (+ (/ (* u4 u4)
200
+					     u3)
201
+					  (/ (* *sos* *sos* u1 u3)
202
+					     (* *kappa* (- *rl* u3)))
203
+					  (/ (* u3 (sigma (phi vars)))
204
+					     *rl*))))))
205
+
206
+;; источниковый член в ячейке
207
+(defmethod source ((vars conservative-variables))
208
+  (let ((d (- ug ul))
209
+	(sg (* -0.5 *rl* d d (signum d))))
210
+    (make-instance 'conservative-variables
211
+		   :var-vector (vector 0 sg 0 -sg))))
212
+
213
+;; градиенты
214
+(defun minmod (a b)
215
+  (if (< (abs a) (abs b)) a b))
216
+
217
+(defun grad-x (grid indx order)
218
+  (let* ((i (x-index indx))
219
+	 (j (y-index indx))
220
+	 (mid-cell (cell grid indx))
221
+	 (left-cell (cell grid (index (1- i) j)))
222
+	 (right-cell (cell grid (index (1+ i) j)))
223
+	 (DL (- (x-coord right-cell) (x-coord mid-cell)))
224
+	 (DR (- (x-coord mid-cell) (x-coord left-cell)))
225
+	 (DUL (- right-cell mid-cell))
226
+	 (DUR (- mid-cell left-cell)))
227
+    (case order
228
+      ((1) (minmod (/ DUL DL)
229
+		   (/ DUR DR)))
230
+      ((2) (/ (+ (/ (* DUL DR)
231
+		    DL)
232
+		 (/ (* DUR DL)
233
+		    DR))
234
+	      (+ DL DR)))
235
+      (otherwise (error "`grad' умеет считать производные лишь до 2го порядка включительно")))))
236
+
237
+(defun grad1-x (grid indx)
238
+  (grad-x grid indx 1))
239
+
240
+(defun grad2-x (grid indx)
241
+  (grad-x grid indx 2))
242
+
243
+;; источниковый член сужения трубки тока за счёт межфазного обмена
244
+(defun ASF (grid indx)
245
+  (let* ((cell (cell grid indx))
246
+	 (phi (phi cell))
247
+	 (rg (rg cell))
248
+	 (ug (ug cell))
249
+	 (ul (ul cell))
250
+	 (C (/ *Cg* (sqrt *kappa*)))
251
+	 (K (/ (* C C rg)
252
+	       *rl*))
253
+	 (N (/ (1- phi)
254
+	       phi))
255
+	 (grads (grad2-x grid indx)))
256
+    (let-conservative grads
257
+      (make-instance 'conservative-variables
258
+		     :var-vector (vector u3
259
+					 (+ (* (- (* C C)
260
+						  (* ug ug))
261
+					       u1)
262
+					    (* 2 ug u2)
263
+					    (* K u3))
264
+					 u4
265
+					 (+ (* C C N u1)
266
+					    (* (- (* K N)
267
+						  (* ul ul))
268
+					       u3)
269
+					    (* 2 ul u4)))))))
270
+
271
+;; выстрелы на границу
272
+(defun shot (grid indx direction)
273
+  (let* ((cell (cell grid indx))
274
+	 (hx (x-size cell))
275
+	 (gx (grad1-x cell indx)))
276
+    (case direction
277
+      ((left) (- cell (* 0.5 hx gx)))
278
+      ((right) (+ cell (* 0.5 hx gx)))
279
+      (otherwise (error 'shot "Неизвестное направление для стрельбы")))))
280
+
281
+(defun left-shot (grid indx)
282
+  (shot grid indx 'left))
283
+
284
+(defun right-shot (grid indx)
285
+  (shot grid indx 'right))
286
+
287
+;; спектральный радиус
288
+(defun radius (cell type)
289
+  (apply #'max (mapcar #'abs (eigen-values-by-cell cell type))))
290
+
291
+;; псевдоскорость
292
+(defun pseudo-vel (grid indx direction)
293
+  (let* ((i (x-index indx))
294
+	 (j (y-index indx))
295
+	 (l-indx (index (1- i) j))
296
+	 (r-indx (index (1+ i) j))
297
+	 (mid-cell (cell grid indx))
298
+	 (left-cell (cell grid l-indx))
299
+	 (right-cell (cell grid r-indx)))
300
+    (case direction
301
+      ((right) (max (radius (right-shot grid indx) 'A)
302
+		    (radius (left-shot grid r-indx) 'A)
303
+		    (radius mid-cell 'Afull)
304
+		    (radius right-cell 'Afull)))
305
+      ((left) (max (radius (right-shot grid l-indx) 'A)
306
+		   (radius (left-shot grid indx) 'A)
307
+		   (radius left-cell 'Afull)
308
+		   (radius right-cell 'Afull))))))
309
+
310
+;(defun border-flux
311
+		     
312
+
313
+
314
+;(defmethod make-step ((grid grid2d) (newgrid grid2d) (tau number))