Browse code

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

А также исправил небольшую багу в определении примитивных переменных через :rg

Dmitrii Kashin authored on 07/05/2015 08:27:48
Showing 8 changed files
... ...
@@ -28,8 +28,10 @@
28 28
 ;;	       (:file "solver2d" :depends-on ("package" "field2d"))
29 29
 ;;	       (:file "march2d" :depends-on ("package" "solver" "grid" "field"))
30 30
 
31
+	       (:file "solver2d-kurganov-tadmor" :depends-on ("grid2d-conservative"))
31 32
 	       (:file "test2d-conservative" :depends-on ("field2d-conservative" "boundary-conditions"))
32 33
 	       (:file "output-functions")
34
+	       
33 35
 	       (:file "march" :depends-on ("grid2d-conservative" "output-functions"))
34 36
 	       ))
35 37
   
... ...
@@ -10,6 +10,7 @@
10 10
 (defparameter *Ro* 287.065476372)
11 11
 (defparameter *T* (/ (* *Cg* *Cg*)
12 12
 		     (* *Ro* *kappa*)))
13
+(defparameter *sigma* 0) ;; добавочное давление
13 14
 
14 15
 ;; конфигурация программы
15 16
 
... ...
@@ -18,3 +19,5 @@
18 18
 
19 19
 (defparameter *number-of-conservative-variables* 4)
20 20
 
21
+(defparameter *newton-max-iter* 50)
22
+(defparameter *newton-est-error* 1.0d-8)
... ...
@@ -78,6 +78,9 @@
78 78
 
79 79
 ;; вспомогательные функции
80 80
 
81
+(defun curryr (function &rest args)
82
+  (lambda (x) (apply function x args)))
83
+
81 84
 (defun filter (&rest args)
82 85
   (apply #'remove-if-not args))
83 86
 
... ...
@@ -1,5 +1,7 @@
1 1
 (in-package :common-lisp)
2 2
 
3
+(load "~/quicklisp/setup.lisp")
4
+
3 5
 (unless *load-truename*
4 6
   (error "This file must be LOADed to set up CL-Frolov."))
5 7
 
... ...
@@ -9,10 +11,12 @@
9 9
 
10 10
 (push (merge-pathnames *clf-home* "systems/") asdf:*central-registry*)
11 11
 (push "/usr/share/common-lisp/source/alexandria/" asdf:*central-registry*)
12
+;(push "/home/freehck/quicklisp/dists/quicklisp/software/gsll-20140211-git" asdf:*central-registry*)
12 13
 
13 14
 (asdf::load-asd "clf.asd")
14 15
 ;(asdf::load-asd "alexandria.asd")
15 16
 
16 17
 (asdf:load-system :alexandria)
18
+(asdf:load-system :gsll)
17 19
 (asdf:load-system :clf)
18 20
 
... ...
@@ -1,5 +1,8 @@
1 1
 (in-package :common-lisp)
2 2
 
3 3
 (defpackage :clf
4
-  (:use :common-lisp :alexandria)
5
-  )
4
+  (:use :common-lisp
5
+	:alexandria
6
+	:gsll)
7
+  ;; следующие символы конфликтовали между alexandria и gsll
8
+  (:shadow mean variance median standard-deviation factorial))
6 9
new file mode 100644
... ...
@@ -0,0 +1,180 @@
0
+(in-package :clf)
1
+
2
+;; определение шага по времени
3
+
4
+(defmethod max-timestep-constant ((grid grid2d-cv) (method (eql :lax-frederichs)) cu)
5
+  5.e-8)
6
+
7
+;; функция-добавок к давлению
8
+(defun sigma (phi)
9
+  phi *sigma*)
10
+
11
+;; производная функции sigma
12
+(defun Dsigma (phi)
13
+  phi 0)
14
+
15
+
16
+(defun eigen-values-by-cell (cell)
17
+  (let* ((phi (phi cell))
18
+	 (rg (rg cell))
19
+	 (ug (ug cell))
20
+	 (ul (ul cell))
21
+	 (rl *rl*)
22
+	 (cg (/ *Cg* (sqrt *kappa*)))
23
+	 (cl (sqrt (/ (+ (* rg cg cg)
24
+			 (* phi (- (sigma phi)
25
+				   rl ul ul)))
26
+		      (* phi rl))))
27
+	 (eigen-values/A ; !!!
28
+	  (list (abs (- ug cg))
29
+		(abs (+ ug cg))
30
+		(abs (+ ul cl))
31
+		(abs (- ul cl))))
32
+	 (eigen-values/sources ; !!!
33
+	  (list (* -1
34
+		   (abs (- ug ul))
35
+		   (+ (/ *rl*
36
+			 (* phi rg))
37
+		      (/ 1
38
+			 (- 1 phi))))))
39
+	 (suggestion (* 5 cg))
40
+	 (left-init (- (min ul (- ug cg) (+ ug cg))
41
+		       suggestion))
42
+	 (right-init (+ (max ul (- ug cg) (+ ug cg))
43
+			suggestion))
44
+	 
45
+	 (X (- (sigma phi)
46
+	       (* (1- phi)
47
+		  (Dsigma phi))))
48
+	 (Z (* (1- phi) rg))
49
+	 (R (+ (* Z cg cg)
50
+		(* X phi)))
51
+	 (c1 (* phi rl))
52
+	 (c2 (* -2 c1 (+ ug ul)))
53
+	 (c3 (- (* c1 (+ (expt (+ ug ul) 2)
54
+			 (* 2 ug ul)
55
+			 (* -1 cg cg)))
56
+		R))
57
+	 (c4 (- (* 2 ug R)
58
+		(* 2 c1 ul
59
+		   (- (* (+ ug ul) ug)
60
+		      (* cg cg)))))
61
+	 (c5 (- (* X cg cg phi)
62
+		(* R ug ug)
63
+		(* c1 ul ul (- (* cg cg)
64
+			       (* ug ug))))))
65
+    (flet ((frolov-equation (x)
66
+	     (- (* phi
67
+		   (- (* rl (expt (- x ul) 2))
68
+		      (- (sigma phi)
69
+			 (* (- 1 phi) (Dsigma phi))))
70
+		   (- (expt (- x ug) 2)
71
+		      (expt cg 2)))
72
+		(* (expt cg 2)
73
+		   (- 1 phi)
74
+		   rg
75
+		   (expt (- x ug) 2))))
76
+	   (frolov-equation-derivative (x)
77
+	     (- (+ (* 2 phi rl
78
+		      (- x ul)
79
+		      (- (expt (- x ug) 2)
80
+			 (expt cg 2)))
81
+		   (* 2 phi
82
+		      (- x ug)
83
+		      (- (* rl (expt (- x ul) 2))
84
+			 (- (sigma phi)
85
+			    (* (- 1 phi) (Dsigma phi))))))
86
+		(* 2 rg
87
+		   (expt cg 2)
88
+		   (- 1 phi)
89
+		   (- x ug))))
90
+	   (frolov-equation-fdf (x)
91
+	     (values (frolov-equation x)
92
+		     (frolov-equation-derivative x)))
93
+	   (find-root (&optional initial (method +newton-fdfsolver+))
94
+	     (let* ((max-iter *newton-max-iter*)
95
+		    (solver (make-one-dimensional-root-solver-fdf
96
+			     method 'frolov-equation 'frolov-equation-derivative 'frolov-equation-fdf initial)))
97
+	       (loop for iter from 0
98
+		  for oldroot = initial then root
99
+		  for root = (progn (iterate solver) (solution solver))
100
+		  while (if (< iter max-iter)
101
+			    (not (root-test-delta root oldroot *newton-est-error* 0.0d0))
102
+			    (error "Метод Ньютона превысил допустимое число итераций."))
103
+		  ;;для отдалки
104
+;		  do (format t "~d~6t~g ~34t~10,6g~&"
105
+;			     iter root (- root oldroot))
106
+		  finally (progn
107
+;			    (format t "~d~6t~g ~34t~10,6g~&~%"
108
+;				    iter root (- root oldroot))
109
+			    (return root))))))
110
+
111
+    (let ((eigen-values/Afull ; !!!
112
+	   (progn
113
+	     ;;(format t "rg: ~d~%" rg)
114
+	     ;;(format t "phi: ~d~%" phi)
115
+	     ;;(format t "ul: ~d~%" ul)
116
+	     ;;(format t "ug: ~d~%" ug)
117
+	     ;;(format t "rl: ~d~%" rl)
118
+	     ;;(format t "cg: ~d~%--~%~%" cg)
119
+	     ;;(format t "X: ~d~%" X)
120
+	     ;;(format t "Z: ~d~%" Z)
121
+	     ;;(format t "R: ~d~%" R)
122
+	     (let* ((root-1 (find-root left-init))
123
+		    (root-2 (find-root right-init))
124
+		    (coefficients (list c1 c2 c3 c4 c5))
125
+		    (quadro-coeffs (divide-polynomial coefficients root-1 root-2)))
126
+	       (multiple-value-bind (root-3 root-4)
127
+		   (apply #'solve-quadratic-complex
128
+			  (mapcar (curryr #'coerce 'double-float) quadro-coeffs))
129
+		 ;;(format t "root: ~d~%" root-1)
130
+		 ;;(format t "root: ~d~%" root-2)
131
+		 ;;(format t "coefficients: ~a~%" coefficients)
132
+		 ;;(format t "quadro-coeffs: ~a~%" quadro-coeffs)
133
+		 ;;(format t "root: ~d~%" root-3)
134
+		 ;;(format t "root: ~d~%" root-4)
135
+		 (list root-1 root-2 root-3 root-4))))))
136
+      
137
+      (values eigen-values/A
138
+	      eigen-values/sources
139
+	      eigen-values/Afull)))))
140
+
141
+
142
+;; чертовски важная штука, без неё reader macro #m не заработает нормально
143
+;; впрочем, я реализовал деление полинома сам, так как в gsll я так и не нашёл нужной функции
144
+(setf grid:*default-grid-type* 'grid:foreign-array)
145
+
146
+;; упрощённая функция деления полинома: производит деление столбиком.
147
+;; coefficients - это список, начинающийся с самой большой степени.
148
+;; деление производится на (x - divisor). Можно указать множество делителей.
149
+;; я её написал потому что не нашёл соответствующую функцию в gcl
150
+(defun divide-polynomial (coefficients &rest divisors)
151
+  (cond ((emptyp divisors) coefficients)
152
+	(t (let ((divisor (first divisors)))
153
+	     (apply #'divide-polynomial
154
+		    (reverse (rest (let ((x (reduce #'(lambda (acc coef)
155
+					       ;;(format t "acc: ~a~%" acc)
156
+					       ;;(format t "coef: ~a~%--~%" coef)
157
+					       (if (emptyp acc)
158
+						   (list coef)
159
+						   (cons (+ coef (* divisor (first acc))) acc)))
160
+					   coefficients
161
+					   :initial-value nil)))
162
+				     ;;(format t "!!! acc: ~a~%" x)
163
+				     x)))
164
+		    (rest divisors))))))
165
+
166
+;; небольшие тесты
167
+
168
+;(defparameter *eigen-cell* (make-instance 'conservative-variables
169
+;					  :phi 0.01
170
+;					  :ug 25
171
+;					  :ul 400
172
+;					  :rg 0.001
173
+;					  ))
174
+;
175
+;(eigen-values-by-cell *eigen-cell*)
176
+						
177
+(defmethod max-timestep/cu ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
178
+  (noop))
179
+
... ...
@@ -3,7 +3,7 @@
3 3
 (defclass conservative-variables ()
4 4
   ((var-vector :initform (make-array *number-of-conservative-variables*)
5 5
 	       :accessor var-vector)
6
-   (headers/conservative-variables :initform '(P phi ug ul)
6
+   (headers/conservative-variables :initform '(P rg phi ug ul)
7 7
 				   :allocation :class)
8 8
    ))
9 9
 
... ...
@@ -13,6 +13,7 @@
13 13
 	(number (setf (var-vector vars) (make-array *number-of-conservative-variables*
14 14
 						    :initial-element number)))
15 15
 	((or phi ug ul P rg)
16
+	 ;(format t "Дано:~%P: ~d~%rg: ~d~%phi: ~d~%ug: ~d~%ul: ~d~%~%" P rg phi ug ul)
16 17
 	 (let ((primitive-vars (make-instance 'primitive-variables
17 18
 					      :P P :rg rg :phi phi :ug ug :ul ul)))
18 19
 	   (setf (var-vector vars)
... ...
@@ -5,7 +5,7 @@
5 5
 (defclass primitive-variables ()
6 6
   ((var-vector :initform (make-array *number-of-primitive-variables*)
7 7
 	       :accessor var-vector)
8
-   (headers/primitive-variables :initform '(P phi ug ul)
8
+   (headers/primitive-variables :initform '(P rg phi ug ul)
9 9
 				:allocation :class)))
10 10
 
11 11
 (defmethod initialize-instance :after ((vars primitive-variables) &key P phi ug ul rg number var-vector)
... ...
@@ -16,7 +16,7 @@
16 16
 	((and phi ug ul (or P rg))
17 17
 	 (progn
18 18
 	   (when P (setf (P vars) P))
19
-	   (when rg (setf (rg vars) (rg->P rg)))
19
+	   (when rg (setf (rg vars) rg))
20 20
 	   (setf (phi vars) phi)
21 21
 	   (setf (ug vars) ug)
22 22
 	   (setf (ul vars) ul)))