Browse code

Метод Курганова-Тадмора: реализовал расчёт шага по времени и собственно сам солвер

Dmitrii Kashin authored on 11/05/2015 01:47:11
Showing 4 changed files
... ...
@@ -21,3 +21,8 @@
21 21
 
22 22
 (defparameter *newton-max-iter* 50)
23 23
 (defparameter *newton-est-error* 1.0d-8)
24
+
25
+(defparameter *multi-root-iter* 1000)
26
+(defparameter *multi-root-est-error* 1.0d-8)
27
+
28
+(defparameter *constant-timestep* .5d-8)
... ...
@@ -262,33 +262,28 @@
262 262
 (defgeneric fixed (grid index border)
263 263
   (:documentation "Граничное условие фиксированных значений"))
264 264
 
265
-
266
-
267
-
268
-;; --------------------------------------------------
269
-
270 265
 (defgeneric border-flux (grid indx tau border)
271
-  (:documentation "Count flux through the border trought the indx's
272
-  cell from grid in time tau. Border can be: :left :right :top
273
-  :down"))
274
-
275
-(defgeneric update-boundary-cells (grid)
276
-  (:documentation "Set boundary cells according to boundary conditions"))
277
-
278
-(defgeneric flux-vector (cell)
279
-  (:documentation "Return flux vector associated with appropriate
280
-  conservative variables"))
266
+  (:documentation "Считает поток через границу ячейки. Граница может
267
+быть одним из значений: :left :right :top :down"))
281 268
 
282
-(defgeneric flux-vector (cv)
283
-  (:documentation "Returns flux vector by conservative-variables vector"))
284
-
285
-(defgeneric max-timestep (grid method)
286
-  (:documentation "Count max timestep for a known method.
287
-Method could be :lax-frederichs"))
269
+(defgeneric max-timestep (grid method cu)
270
+  (:documentation "Считает максимальный шаг по времени для выбранного
271
+метода `method'. Метод может быть одним из следующих
272
+значений: :lax-frederichs :kurganov-tadmor"))
288 273
 
289 274
 (defgeneric make-step (grid newgrid tau)
290
-  (:documentation "Take grid and make a time-step tau. The result is
291
-  written into newgrid"))
275
+  (:documentation "Делает в сетке `grid' шаг по времени величиной
276
+  `tau'. Результат записывается в `newgrid'"))
292 277
 
278
+(defgeneric flux (vars)
279
+  (:documentation "Возвращает потоковый вектор, рассчитанный по
280
+  переменным `vars'"))
293 281
 
282
+(defgeneric source (vars)
283
+  (:documentation "Возвращает вектор источниковых членов, рассчитанный
284
+  по переменным `vars'"))
294 285
 
286
+;; --------------------------------------------------
287
+
288
+(defgeneric update-boundary-cells (grid)
289
+  (:documentation "Set boundary cells according to boundary conditions"))
... ...
@@ -2,8 +2,8 @@
2 2
 
3 3
 (defun march (grid &key output-every-step output-every-time finish-time finish-step
4 4
 		     print-out-border suppress-export
5
-		     (max-timestep #'max-timestep-constant)
6
-		     (make-step #'make-step/wo-sources)
5
+		     (max-timestep #'max-timestep)
6
+		     (make-step #'make-step)
7 7
 		     (cu nil))
8 8
   ;; Проверка начальных условий и вывод информации о процессе
9 9
   (format t "Маршевая функция запущена.~%")
... ...
@@ -149,15 +149,6 @@
149 149
 ;;(eigen-values-by-cell *eigen-cell* 'S)
150 150
 ;;(eigen-values-by-cell *eigen-cell* 'Afull)
151 151
 
152
-;; определение шага по времени
153
-
154
-(defmethod max-timestep-constant ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
155
-  5.e-8)
156
-
157
-;; к сожалени, пока нет описания, как это делать
158
-(defmethod max-timestep/cu ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
159
-  (noop))
160
-
161 152
 ;; потоковый член в ячейке.
162 153
 (defmethod flux ((vars conservative-variables))
163 154
   (let-conservative vars
... ...
@@ -279,7 +270,108 @@
279 279
 		   (radius left-cell 'Afull)
280 280
 		   (radius right-cell 'Afull))))))
281 281
 
282
-		     
283 282
 
283
+;; определение шага по времени
284
+
285
+(defmethod max-timestep ((grid grid2d-cv) (method (eql :kurganov-tadmor)) cu)
286
+  (if (not cu)
287
+      *constant-timestep*
288
+      (let ((tau 1000))
289
+	(loop for i to (cells-max-x grid) do
290
+	     (loop for j to (cells-max-y grid) do
291
+		  (let* ((index (index i j))
292
+			 (cell (cell grid index))
293
+			 (h (max (x-size cell)
294
+				 (y-size cell)))
295
+			 (a+ (pseudo-vel grid index 'right))
296
+			 (a- (pseudo-vel grid index 'left))
297
+			 (new-tau (/ (* h cu) (max a+ a-))))
298
+		    (when (< new-tau tau)
299
+		      (setf tau new-tau)))))
300
+	tau)))
301
+
302
+;; потоки через границы ячейки
303
+
304
+(defmethod border-flux ((grid grid2d-cv) (indx+ cell-index) (tau number) (border (eql :left)))
305
+  (let* ((indx- (index (1- (x-index indx+))
306
+		       (y-index indx+)))
307
+	 (u+ (left-shot grid indx+))
308
+	 (u- (right-shot grid indx-))
309
+	 (F- (flux u+))
310
+	 (F+ (flux u-))
311
+	 (a (pseudo-vel grid indx- 'right)))
312
+    (cvar* 0.5
313
+	   (cvar- (cvar+ F+ F-))
314
+	   (cvar* a (cvar- u+ u-)))))
315
+
316
+(defmethod border-flux ((grid grid2d-cv) (indx cell-index) (tau number) (border (eql :right)))
317
+  (let ((indx (index (1+ (x-index indx))
318
+		     (y-index indx))))
319
+    (border-flux grid indx tau :left)))
320
+
321
+;; солвер
284 322
 
285
-;(defmethod make-step ((grid grid2d) (newgrid grid2d) (tau number))
323
+(defmethod make-step ((grid grid2d-cv) (newgrid grid2d-cv) (tau number))
324
+  (labels
325
+      ((predictor (grid tau)
326
+	 (let ((newgrid (clone grid)))
327
+	   (loop for j to (cells-max-y grid) do
328
+		(loop for i to (cells-max-x grid)
329
+		   for dF+ = (border-flux grid (index i j) tau :right)
330
+		   for dF- = (border-flux grid (index i j) tau :left) then dF+
331
+		   do
332
+		     (let* ((dF (cvar- dF+ dF-))
333
+			    (u (cell grid (index i j)))
334
+			    (hx (x-size u))
335
+			    (s (source u)))
336
+		       (flet ((eq-system-non-linear (x1 x2 x3 x4)
337
+				(let* ((x (make-instance 'conservative-variables
338
+							 :var-vector (vector x1 x2 x3 x4)))
339
+				       (result (cvar- (cvar+ (cvar/ (cvar- x u)
340
+								    tau)
341
+							     (/ dF hx))
342
+						      (* 0.5 (cvar+ s (source x))))))
343
+				  (let-conservative result
344
+				    (values u1 u2 u3 u4)))))
345
+			 (let-conservative u
346
+			   (let ((max-iter *multi-root-iter*)
347
+				 (solver (make-multi-dimensional-root-solver-f 
348
+					  method #'eq-system-non-linear
349
+					  (grid:make-foreign-array
350
+					   'double-float :initial-contents (list u1 u2 u3 u4)))))
351
+			     (loop for iter from 0
352
+				;;for fnval = nil then (function-value solver)
353
+				;;for argval = nil then (solution solver)
354
+				while (if (< iter max-iter)
355
+					  (or (zerop iter)
356
+					      (not (multiroot-test-residual solver *multi-root-est-error*))))
357
+				do (iterate solver)
358
+				finally (setf (cell newgrid (index i j))
359
+					      (make-instance 'conservative-variables
360
+							     :var-vector (solution solver))))))))))
361
+	   newgrid))
362
+       (corrector (grid predgrid tau)
363
+	 (let ((midgrid (clone grid))
364
+	       (newgrid (clone grid)))
365
+	   (loop for j to (cells-max-y grid) do
366
+		(loop for i to (cells-max-x grid) do
367
+		     (setf (cell midgrid (index i j))
368
+			   (cvar* 0.5 (cvar+ (cell grid (index i j))
369
+					     (cell predgrid (index i j)))))))
370
+	   (loop for j to (cells-max-y grid) do
371
+		(loop for i to (cells-max-x grid)
372
+		   for index = (index i j)
373
+		   for dF+ = (border-flux midgrid index tau :right)
374
+		   for dF- = (border-flux midgrid index tau :left) then dF+
375
+		   do
376
+		     (let ((S (cvar* 0.5 (cvar+ (source grid index)
377
+						(source predgrid index))))
378
+			   (u (cell grid (index i j)))
379
+			   (hx (x-size u))
380
+			   (dF (cvar- dF+ dF-)))
381
+		       (setf (cell newgrid index)
382
+			     (cvar- (cvar+ u (cvar* tau S))
383
+				    (cvar/ (cvar* tau dF)
384
+					   hx))))))
385
+	   newgrid)))
386
+    (corrector grid (predictor grid tau) tau)))